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Abstract 


Improved measurements of prectpitatton wi„ aid our undoing ofdte ro ,e ofiatent 

heattng on global emulations. Spacebonte meteorological senses such as the pianned 

P on radar and microwave radiometers on the Tropical Rainfall Measurement 

Mtsston (TRMM) provide for the firs, „me a comprehensive means of making these global 

measurements. Pre-TRMM activities tnclude development ofptecipitafion algorithms 

ustng ex, sung satelltte dam, computer simulations, and measurements ftom limited aircraft 

campaigns. Since the TRMM radar will be the fts, spacebome ptecipitation radar, them is 

limited expertence with such measurements, and only recendy have atrbotne radars 

become avatlable that can attempt to address the issue of the limitations of a spacebome 

radar. Tftere are many qU es„ons regarding how much attenuation occurs in vanous cloud 

types and the effect of cloud vertical motions on the estimation of precipitation rates. The 

program betng developed by NASA GSFC will provide dah, useful for testing both 

ram-tetneval algortthms and the impoftance of venical motions on the nun measurements 

TTa Purpose of this report is to describe the design and development of real-time 

embedded parallel algortthms used by EDOP to extract reflectivtty and Doppler products 

(velocity, spectrum w.dth, and signal-, o-notse ratio, as the firs, step in the aforemen, toned 
goals. 
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Chapter 1 - Introduction 


1.0 Introduction 

Global climate models depend upon knowledge of the atmospheric variables on a 
worldwide scale. Though satellites have greatly improved the inputs to these models, 
much remains to be done before inputs are adequate for long-range forecasts and study of 
major global changes. In particular, we can only crudely model the major heat exchange 
between atmosphere and surface that takes place in tropical rainfall (Simpson, 1988). 
Current models for global circulation depend largely on sparsely-spaced upper-air 
soundings and surface measurements at infrequent intervals in time and space. Moreover, 
lack of most meteorological parameters on model grid scales (200 km) at the surface and 
aloft over the oceans represents a major gap in knowledge. Clearly this information is a 
requirement to improve global circulation models. 


Tropical rainfall is critical to the distribution of life-giving water throughout our Earth 
system. Over two-thirds of the worldwide precipitation falls in the Tropics, releasing 
latent heat and thereby powering the global atmospheric circulation that shapes our world 
climate. English scientist G. Hadley (18th century) was the first to understand the central 
role played by tropical rainfall in shaping Earth’s climate. In brief, cool low-level moist air 
flowing toward the equator replaces heated equatorial air as it ascends and moves toward 
the poles at high altitudes. Stored solar energy is released as latent heat when the water 
vapor condenses into cloud clusters or Hadley cells reaching hundreds of kilometers in 
horizontal extent. Of these, perhaps one in ten deepens into a tropical cyclone. The solar 
energy released as latent heat by raining tropical clouds is the main source of energy for 
global atmospheric motions. In fact, low-level Hadley circulation carries evaporated water 
from over one-half the earth's surface into the Inter Tropical Convergence Zone (ITCZ) 
where it is precipitated by raining tropical systems. 
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Rainfall in the Tropics additionally plays a fundamental role in the El Nino climate 
anomalies whose disruptions to worldwide precipitation patterns trigger floods and 
drought around the globe. It’s four- to five-year-average cycle is associated with 
anomalous warming of equatorial surface waters extending from the Eastern Pacific to the 
coast of South America. This warming is nearly always accompanied by large-scale shifts 
in atmospheric pressure that go from the East to the West Pacific and back in an 
oscillatory manner often referred to as the Southern Oscillation. The El Nino and 

Southern Oscillation are merely two manifestations of the same climatic disturbance— the 
ENSO event. 


The importance of determining, at least climatologically, the rainfall in the Tropics is the 
subject of many papers and reports leading to the Tropical Rainfall Measurement Mission 
(TRMM). Lau et al. state that latent heating by precipitation causes about 40% of the 
energy flux to and from the atmosphere. This takes place largely in poorly instrumented 
tropical regions: oceans with little ship traffic and less-developed land areas with few 
weather stations. The information on this rainfall is so sparse that we do not even have 
adequate climatological data, much less synoptic data. The purpose of TRMM is to 
provide a first look at this climatology globally below 35° latitude. 

Whereas TRMM will be able to measure rainfall rates on a climatological basis, it cannot 
provide wind velocities over the same tropical environment. Observations of horizontal 
winds and convection on a global scale are necessary for advancing our knowledge of 
large-scale atmospheric circulation and climate dynamics, thereby improving our 
understanding of the global biogeochemical and hydrologic cycles. Wind measurements 
are particularly important in the Tropics where there are few observing stations— indeed, 
it has been suggested (NASA, 1987) that wind profiles are the single most important new 
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data source because of their potential for dramatically advancing our skill at Numerical 
Weather Prediction (NWP). 

1.1 Earth System Interactions 

The atmosphere, oceans, and land biomass systems are closely coupled through the 
unifying role of global hydrology. Because of this, studies of rainfall are essential to gain 
deeper understanding of our Earth system. 

General Circulation Models (GCM) simulate general atmospheric features to investigate 
the climate variations that take place on a variety of temporal scales. Existing GCM 
predictions are extremely sensitive to the assumed vertical distribution of latent heating by 
equatorial cloud systems; however, the vertical distribution of cloud layers is poorly 
observed. This casts doubts on our understanding of the role played by clouds in the 
climate system and their effects on the surface-radiation balance and profile of heat 
absorption and re-radiation throughout the atmosphere. 

The response of cloud systems to environment is an important link in a chain of processes 
responsible for monsoons, ENSO events, and other climate variations. Numerical models 
of precipitation systems provide essential insights to the interaction of clouds with each 
other, their surroundings, and land/ocean surfaces. They are needed to convert rain-rates 
derived utilizing radiometric temperatures from spacebome microwave and radiometers 
into latent-heating profiles. Such models further aid computations of the total energy 
budget, which consists of radiative and latent heating effects. Though poorly observed, 
cloud profiles strongly impact the hydrologic and carbon cycle. Clearly, improved 

observations of clouds are critical for updating global circulation models used in climate 
prediction. 
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Ocean-atmosphere interactions likewise render a key cole in the Earth system. For 
example, as horizontal winds stress the sea surface, oceanic circulations develop 
influencing the transport of heat across the ocean-atmosphere boundary. Fresh water 
from precipitation further alters sea-surface temperatures as do clouds that shield the 
ocean surface from solar radiation. The result is a complex web of interaction between the 
wind and water that acts to regulate our climate. 


In equatorial forest regions, tropical rainfall directly influences the terrestrial water and 
energy balance. Widespread cloudiness shields tropical forests from direct sunlight 
encouraging high growth rates and a large standing biomass. Most of the water received 
as rainfall is returned to the atmosphere in ten-day cycles through evaporation and 
transpiration, leading us to view the tropical forest and atmosphere as a single water and 
energy regulation system. Forests and soils are also a major source of many atmospheric 
trace constituents (e.g., carbon monoxide, methane, and the hydroxyl radical) that are 
taken into the lower atmosphere during convective rainstorms. This exchange provides a 
direct link between rainfall and the global biogeochemical cycles that regulate life on 
Earth. 

1.2 What now? 

Existing spacebome atmospheric sensors are passive in the visible, infrared (IR), and 
microwave parts of the spectrum, and the measured radiometric temperatures are 
vertically integrated quantities. The visible and IR instruments sense upwelling radiation 
near cloud tops. Passive-microwave instruments provide integrated liquid water and 
water vapor estimates. Thus, interpretations of these temperatures to retrieve physical 
quantities are often difficult and frequently necessitate the use of complementary 
instruments to provide vertical structure because of the inherent lack of passive, vertical- 
profiling capability. Now we can look forward to development and deployment of 
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spacebome radars to provide vertical precipitation structure for meteorology. The first of 
these will be on the TRMM, scheduled for launch in 1997, and will provide vertical 
precipitation profiles using the radar reflectivity. Such measurements can then be 
combined with cloud models to estimate profiles of latent heating. 

The basic problem with existing studies is that they have been tested with limited aircraft 
data. Intensive aircraft campaigns are necessary to validate cloud models and algorithms 
as well as provide case studies and reference measurements for their spacebome 
counterparts. To test them, both radars and radiometers must fly above thunderstorms 
and precipitation systems in regions where surface measurements of precipitation are 
available. Whereas surface measurements provide rain rates and total rain amounts 
directly, aircraft or satellite measurements provide only indirect indications of the rain rate 
from the backscattered power (or radar reflectivity). These estimates are based on 
empirical relations and can be ambiguous. For example, they assume that vertical air 
motions in clouds are negligible, when in fact they can be substantial. Existing rain- 
retrieval algorithms fail to utilize knowledge of vertical wind motions because they are 
unavailable. While it has been argued that these errors are minimized in the monthly- 
averaged data such as will be provided by TRMM, it is desirable to more closely examine 

the effect of vertical motions on rain-rate estimates using aircraft instrumented with radar 
and radiometers. 

An X-band Doppler radar, the EDOP developed by the National Aeronautics and Space 
Administration's Goddard Space Flight Center (NASA GSFC), has now begun flying at 
high altitude on the ER-2 aircraft (Heymsfield et al„ 1989; Heymsfield et al„ 1993). This 
system will make possible the testing of algorithms for rain retrieval, the concepts for 
radar wind measurement, and the validation of existing cloud models. 
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1.3 Purpose of this paper 

The purpose of this Masters project is to document the design and development of real- 
time embedded parallel-processing algorithms executed by the EDOP Data System. These 
algorithms reduce raw radar measurements to reflectivity and Doppler products (velocity, 
spectrum width, and signal-to-noise ratio) resulting an effective compression of the raw 
data exceeding 400: 1 . 

To provide the reader with necessary background material. Chapter Two introduces 

EDOP's scientific objectives and compares the role of EDOP to existing airborne radars in 
meeting these needs. 

EDOP is understandably a complicated system and discussion of its signal processing 
algorithms is incomplete without an overview of the system design. Chapter Three 
examines the RF and digital hardware in limited detail to give the reader an appreciation of 
the interaction between various system components during operation and the manner in 
which they constrain algorithm development. 

The design and development of the reflectivity and Doppler signal-processing algorithms 
are investigated in detail in Chapter Four. This chapter forms the basis for this report and 
additionally derives estimates of the measurement accuracies. 


Chapter Five concludes this report by exploring the remaining tasks needing completion 
before EDOP fully meets its specified operational state. Two remaining tasks of particular 
importance include nadir antenna stabilization and development of an automatic gain 

control to prevent receiver saturation due to the large dynamic range of received weather 
echoes. 
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Chapter Two - Background 


2.0 Historical Digression 

The earliest known treatises on meteorology date back to Aristotle’s (384-322 B.C.) 
publication of Meteorologica. Though his pupil, Theophrastus, continued the work of 
writing about winds and weather signs, progress was excruciatingly slow for the next 2000 
years largely due to a lack of instruments for observing the primary measurables. The 
early seventeenth century heralded rapid growth in meteorological research due to the 
invention of the thermometer (Galileo, 1607) and barometer (Evangelista Torricelli, 1643) 
as well as the discovery of Boyle’s Law (1659). Less than forty years later, Edmund 
Halley (1696) ushered the era of Earth System Science through his attempt to explain 
general circulation by variable surface heating. Further attempts were made to explain 
atmospheric motions during the 1700's (Hadley, 1735; Lavoisier, 1783), but it wasn’t until 
John Dalton realized in 1800 the relation between expanding air and atmospheric 
condensation that the physical basis of modem meteorology was established. Between 
1800 and 1815, J. B. Lamarck, with the help of Lavoisier, Laplace, and others, created the 
first international compilation of weather observations. These synoptic observations 
affirmed the existence of characteristic patterns of pressure, temperature, winds, and, 
moreover, empirical rules for their development, movement, and the consequent weather 
changes. Such knowledge was of immediate benefit to the large number of sailing ship 
captains who invested in M. F. Maury's 1848 publication of global wind-field maps. The 
development of rubber balloons (1890's) led to the frequent sampling of the upper 
atmosphere by France, Germany, and England. Observations of winds aloft could then be 
made by watching the motion of the balloons, provided they didn’t enter the clouds! 

The first verified report of a precipitation-radar echo was by a 10 cm system that tracked a 
shower to a distance of seven miles off the English Coast on 20 February, 1941. Some 
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very weak echoes were thought to be attributed to storms as early as 1938; however, the 
detection of hydrometeors was severely limited by the fact that available methods of RF 
power generation resulted maximum frequencies between 200 and 400 MHz. The 
subsequent development of the highly secret British magnetron expanded the boundary to 

3000 MHz (S-band) and later 10,000 MHz (X-band) inspiring the genesis of radar 
meteorology. 

2.1 Scientific Objectives 

The EDOP radar's major justification is to provide a means to study the dynamic and 
hydrometeor structure of thunderstorms and mesoscale convective systems using a 
comprehensive instrument package on the ER-2. The aim is to improve the 
understanding of thunderstorms and larger scale or mesoscale convective systems using 
present and future satellites (e.g. GOES, TRMM, EOS-MIMR). Specifically, EDOP's 
scientific objectives are to provide radar data in conjunction with the multi-frequency 
measurements from other ER-2 complementary instruments to: 

• Better understand the structure of deep isolated thunderstorms, 

• Determine the structure and dynamics of both the overshooting and anvil 
portions of individual thunderstorms, 

• Define the mesoscale structure and dynamics of mesoscale convective systems 
(MCS) including squall lines and mesoscale convective complexes, 

• Better understand the cloud microphysics in isolated thunderstorms and 
MCS's, 

• Test proposed satellite passive- and active-microwave precipitation-estimation 
algorithms. 
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2.1.1 Structure of deep isolated thunderstorms 

It is of interest to utilize satellite IR imagery for nowcasting of thunderstorms, 
precipitation estimation, and studies of thunderstorm anvil structure. Intense midwest 
thunderstorms frequently exhibit peak vertical-wind motions of 40-50 m/s and 
reflectivities exceeding 60 dBZ with large hail in the vicinity of the updraft core. 
Numerous studies have probed the updraft structure utilizing multiple ground-based 
Doppler radars by integration of the mass-continuity equation which uses horizontal wind 
measurements to deduce vertical wind speeds. Yet these vertical-wind speeds estimates 
often have accuracies only to within a few meters per second. Contamination of radial 
velocities by ground clutter and weak signal-to-noise ratios (SNR) result in inadequate 
estimates of mass divergence in the boundary layer leading to poor estimates of vertical 
velocities (w). Interpolation of radial velocities from different radars to a common 3- 
dimensional grid, grid filtering, and estimation of top and bottom boundary conditions 
required for w calculations further degrade vertical velocity accuracy. Moreover, typical 
radar sampling resolutions limit resolvable scales of motion to greater than 5 km 
(Carbone, 1985). A vertically pointing airborne radar can mitigate these problems. 

2.1.2 Anvil structure and dynamics 

Satellite infrared observations of severe thunderstorms frequently show a thermal couplet 
near cloud tops and a "V" shaped region of cold temperatures on the anvil scale 
(Heymsfield et al, 1983). Explanations of these phenomena are controversial; however, it 
is generally accepted that cloud air is warmed as overshooting and negatively buoyant 
updraft air subsides downwind of the cloud top (Heymsfield, 1983; Schlesinger, 1984; 

Adler and Mack, 1986). An alternative explanation (Heymsfield and Blackmer, 1988), 
schematically illustrated in Figure 2-1, is based on the generation of lee waves by 
stratospheric flow over the penetrating cloud tops; however, no insitu data exist to 
confirm the various hypotheses. This inadequate knowledge of vertical-motion structure 
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and precipitation profiles hampers effective utilization of GOES imagery to deduce the 
mixing of penetrating updraft air with stratospheric air. 

Both airborne and ground-based radars typically determine vertical velocity by integrating 
horizontal winds using the mass-continuity equation. The limits of integration require 
knowledge of the top and bottom boundary conditions; however, existing radars regularly 
fail to see the small ice crystals in the tops of storms, and the resulting ambiguous 
estimates of the top boundary introduce error into the vertical-velocity estimates. An ER- 
2 radar can measure the weak cloud top echoes with greater accuracy because of its 
close proximity — often the cloud tops are only 3-5 km below the aircraft. 

In addition to the overshooting portion of the storm, the surrounding anvil is of great 
interest. Anvils produced by deep isolated thunderstorms are frequently extensive — much 
larger than the core region itself. The horizontal extent, depth, and hydrology of the anvil 
are related to storm factors such as updraft intensity, longevity, precipitation efficiency, 
and vertical shear; yet, convection can perturb the mesoscale environment in such a way 
as to feed on the convection itself. To better understand the complete dynamics of deep 
convection, it is necessary to obtain improved observations of air motions within the anvil 
region. Ground-based radars inherently obtain poor samples of the anvil region because of 
its large horizontal extent (hundreds of kilometers). Furthermore, such radars are rarely 
situated directly under the thunderstorm for ideal viewing; instead, they must observe 

distant events. The mobility of an ER-2 radar makes it ideally suited for anvil-region 
studies. 
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2.1.3 Structure of mesoscale convective systems 

Within the past decade, considerable attention has been devoted to probing the two- and 
three-dimensional structure of mesoscale convective systems (Figure 2-2). These systems, 
extending a few hundred kilometers, are often composed of both convective and stratiform 
precipitating regions. Updrafts of 20 m/s or more are common in the convective region 
whereas the weaker stratiform circulations are driven largely by diabatic heating and 
cooling processes such as melting of ice or evaporation of rain. While multiple Doppler 
networks of ground based radars can be used to observe these mesoscale systems, 
airborne platforms are better suited due to their greater mobility. Conventional techniques 
require integration of the mass-continuity equation thereby resulting estimates of vertical 
motion with sometimes larger than desirable uncertainties. Recently, the scanning pattern 
of the NOAA P-3 radar has been changed to utilize an approach similar to that of 
ELDORA. This has provided effective vertical incidence data in limited regions flyable by 
the P-3 (Marks and Houze, 1987). An ER-2 based radar would have the mobility to 
provide similar measurements from directly above convective storms. 


2.1.4 Cloud microphysics 

The EDOP radar should additionally have the capability to deduce the precipitation phase 
of hydrometeors by making polarimetric measurements. The linear depolarization ratio 
(LDR) is conventionally defined as the ratio of received power at two orthogonal 
polarizations from a linearly-polarized transmitted pulse: 

LDR = lOlog— 

Zhh 


Hailstones tend to show no preference to polarization because they tumble during descent 
due to their large size and typically exhibit LDR ratios of -10 dB. Raindrops, however, 
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tend to flatten while falling resulting in ratios less than -30 dB. Snow takes on 
intermediate values. The design of EDOP's antennas and radome provide moderate 
isolation (-28 dB) between the polarizations; nevertheless, this should be sufficient to 
discriminate hail from other hydrometeor types thereby providing additional justification 
for assuming a particular hydrometeor f allspeed. 

2.1.5 Testing of satellite precipitation algorithms 

Ground-based radars offer excellent spatial and temporal coverage over limited regions, 
but often weather systems are beyond range of these radars (e.g. no data exist over the 
oceans). These ground-based measurements consequently provide insufficient 
precipitation estimates for use in global-scale weather prediction, climate, and radiative- 
transfer models. Spacebome meteorological sensors such as the planned precipitation 
radar and microwave radiometers on the Tropical Rainfall Measurement Mission (TRMM) 
provide for the first time a comprehensive means of making these global measurements, 
but much needs to be learned about algorithms combining passive-microwave and radar 
measurements. Necessary pre-TRMM activities include development of precipitation 
algorithms using existing satellite data, computer simulations, and measurements from 
limited aircraft campaigns. Since the TRMM radar will be the first spacebome 
precipitation radar, there is limited experience with such measurements, and only recently 
have airborne radars become available that can attempt to address the issue of the 
limitations of a spacebome radar. Although the operating frequency of EDOP is different 
from that of TRMM, the less-attenuated wavelength should be valuable in this "combined- 
algorithm development and provide necessary ground truth for calibrating spacebome 
radars. The dual-beam approach used by EDOP mil additionally provide a unique 
opportunity to test the dual-beam rainfall-estimation algorithm proposed by Testud 
(1987) for use with airborne and spacebome systems. This method uses a variational 
approach to obtain path attenuation and relates it to the rain rate using a power law. 
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2.2 Existing and proposed airborne radars 

There presently exists an abundance of airborne radars which can be grouped into two 
general categories: precipitation and cloud radars. Precipitation radars are used to 
measure the structure and evolution of actively precipitating and possibly convective 
systems whereas cloud radars penetrate and observe the vertical distribution of cirrus, 
stratocumulus, shallow cumulus, and other radiatively important cloud types. Often these 
systems are designed for specific applications so it is important to confer the uniqueness of 
EDOP in relation to existing and proposed systems (Table 2-1) and show its role in 
advancing our understanding of meteorology. 


The NOAA P-3 radar is mounted on the tail of the P-3 aircraft and scans in a helical 
fashion to efficiently radiate the space around the aircraft. This radar, in operation since 
1983, can sample reflectivity fields and provide both single- and dual-Doppler 
observations from a maximum altitude of about 8 km. In the dual-Doppler mode, the 
aircraft must fly a checkerboard flight track such that alternating flight tracks are normal 
to each other. If the observed storms undergo significant changes during the data 
collection period, the time difference between orthogonal tracks can result in large errors 
in the calculated wind fields (Hildebrand and Moore, 1992). 

The ELectra DOppler RAdar (ELDORA) is an X-band, scanning, dual-beam, stepped- 
frequency system developed jointly by the National Center for Atmospheric Research 
(NCAR) and the French Centre de Recherche en Physique de l’Environnement (CRPE). 
Under an agreement signed in 1990, NCAR supplied the transmitter, receiver, signal 

processor, data system, and aircraft whereas CRPE provided the antennas (Hildebrand, 
Testud, 1992) 
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ELDORA/ASTRAEA Hies between 8 and 10 km on the NCAR Electra making radial 
dual-Doppler velocity and reflectivity measurements of storm structure and kinematics 
with sufficient accuracy to derive thermodynamic properties of the storm such as 
buoyancy, heat and water transport. 

NASA s Jet Propulsion Laboratory began development of the Airborne Rain Mapping 
Radar (ARMAR) in 1986. This DC-8 based, dual-wavelength, polarization diverse, pulse- 
compression system will largely be used to validate TRMM (Im et al„ 1992). It's short 
term goals are to verify the TRMM radar system design, allow development and testing of 
spacebome rain-retrieval algorithms, and serve as a reference for the validation of TRMM 
data after launch in 1997. This dedicated research tool will further serve as a testbed for 
developing future spacebome techniques and technologies for spacebome rain radars. 


An airborne rain-scatterometer/radiometer has been developed by the Radio Research 
Laboratories (CRL) in Japan (Okamoto et al., 1982). This dual-wavelength, X- and K a - 
band system has been operated on a variety of aircraft (e.g. Cessna 404, NASA T-39, 
NASA DC-8, NASA P-3) at altitudes up to 12 km. This system was originally designed 
to study the bright-band structure of rain but has since been used to examine hydrometeor 
growth processes. It should be noted that this is a non-Doppler system. 

The NASA Goddard Space Flight Center is presently developing a W-band (94-GHz), 
polarization-diverse, scanning, Doppler radar to be flown with a host of ER-2 based 
instruments including EDOP, CLS, AMPR, and MIR (Racette, Heymsfield, 1993). The 
Doppler observations will provide insight into the dynamics of clouds. It will further serve 
as a useful test-bed to determine the utility of a future spacebome cloud radar. 


20 


The Universities of Massachusetts and Wyoming have collaborated to produce a fully 
polarimetric, W-band, King-Air based cloud radar (Vali, 1993). This system can be 
pointed either horizontally or vertically perpendicular to the flight path to determine co- 
and cross-polarization power and phase and Doppler velocity. Preliminary test flights 
revealed an unexpected sharp peak in LDR at the melting-layer boundary. It is believed 
that this is the result of an asymmetry in the melt-water distributions on the ice crystals. 

Clearly no one instrument can do it all. A comprehensive package of passive- and active 
microwave sensors such as those on the ER-2 in conjunction with EDOP is ideal for 
study of these complicated phenomena. 


Table 2-1: Existing and Planned Airborne Meteorological Radars 


Instrument Frequency Band Operator 


NOAA 


NOAA-P3 

X 

ELDORA/ 

ASTRAEA 

X 

ARMAR 

K.K,, 


KING-AIR W 


Comments 


scanning, dual- 
Doppler 




NASA/GSFC 


Universities of 
Massachusetts/ 
Wyomin 


polarization- 
diverse, pulse- 
compression, 
TRMM simulator 


| scanning, 
polarization- 
diverse, cloud 
radar, still in 
development phase 


fully polarimetric 
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Chapter 3 - System Design 


3.0 System Overview 

EDOP is an X-band, dual-beam, polarization-diverse weather radar developed for the 
NASA ER-2 aircraft (Heymsfield et. al„ 1989; Heymsfield et. al., 1991). Figure 3-1 
illustrates a typical configuration of EDOP and complementary instruments in the ER-2 
payload bays. The system flies at a nominal altitude of 20 km (70,000 ft.) to map out 
time-height sections of hydrometeor reflectivity as well as vertical and along-track wind 
velocities (Figure 3-2). The nadir antenna measures co-polarization reflectivity, Doppler 
velocity, and spectrum width. Its velocity measurements are used to provide a direct 
indication of vertical air motion once the hydrometeor fall speed is removed. A second 
forward-pointing antenna (30° from nadir) measures cross-polarization reflectivity in 
addition to the previous parameters and can be used to determine the linear depolarization 
ratio (LDR) useful for assessing hydrometeor type and providing additional justification 
for the assumed fall speed. 

Real-time reflectivity and Doppler processing of seven (four linear and three logarithmic) 
channels of data are accomplished by a data system built in-house at Goddard. Two four- 
channel A/D VME 9U cards sample the analog signals at 4 MHz to provide 37.5 m 
resolution in range. Three Martin Marietta LUA-200 linear-mapped array signal- 
processing boards each containing eight AT&T WEDSP32C digital-signal-processors 
provide peak computational speeds of 600 million floating-point operations per second 
(MFLOPs) when clocked at 50 MHz. A Radar Interface Card (RIB) establishes direct 
control of transmitter and receiver functions. Finally, a 68030 VME-based Eltec E-6 
Single Board Computer is used as a host. 
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Figure 3-1 : 


Typical muli t-instrument configuration of ER-2. 



EDOP Measurement Concept 



Figure 3-2: 


EDOP measurement concept (Heymsfield, 1993) 
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3.1 RF Section 

Performance specifications for EDOP and a system block diagram are provided in Table 3- 
1 and Figure 3-3, respectively. In addition, Figure 3-4 illustrates the layout of system 
components in the ER-2 nose cone. 


Table 3-1: EDOP Specifications 


Frequency 

9.72 GHz 

Peak Power 

20 kW 

Duty Cycle 

1% max 

Pulse Length 

0.25, 1.00 us 

PRF 

2200, 4400 Hz 

Antenna Diameter 

0.76 m 

Antenna Beamwidth 

2.9° ] 

First Sidelobe Level 

-30 dB 

Cross-polarization Level 

-38 dB “l 

Receiver Dynamic Ranee 

1 10 dB 

Number of Doppler Channels 

2 

Number of Log Reflectivity Channels 

3 

Nadir-Beam: 


Transmit Polarization 

Horizontal 

Received Polarization 

Co-polarized 

Forward-Beam: 


Transmit Polarization 

Vertical 

Received Polarization 



Co-pol. and 
Cross-pol. 


A 9.66 GHz coherent phase-locked oscillator is used to generate the transmit carrier 
frequency of the system. This RF energy is mixed with the 60 MHz local-oscillator and 
gated by a set of switched PIN diodes before being amplified by a high-gain 20 kW Litton 
air-cooled Traveling Wave Tube (TWT) amplifier and routed through circulators to dual 
high-gain (36 dBi) offset-fed parabolic antennas. The existing mode of operation allows 
power to be radiated simultaneously through both antennas; however, some consideration 
has been given to a future second-mode where power is radiated in alternate pulses 
between the two antennas. 
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A low-noise RF gallium-arsenide pre-amplifier positioned near the antennas limits the 
receiver noise figure to approximately 1 .8 dB. Nadir co-polar and forward co- and cross- 
polarized reflectivity channels pass through logarithmic detectors at the 60 MHz IF and 
then to the data system for real-time processing. Meanwhile, linear synchronous detectors 
beat dual-Doppler RF channels to baseband providing the data system with four (4) I and 
Q channels for Doppler-velocity extraction. Digitally controlled RF and IF attenuators are 
utilized to manipulate receiver gain, but an automatic gain-control (AGC) algorithm has 
not yet been implemented. 

3.2 Data System 

The major objective of real-time processing is the reduction of the raw data to levels that 
can be stored using cost-effective technology. It would be preferable to store raw data 
and perform off-line processing with a variety of algorithms; however, a simple calculation 
illustrates that this goal is not currently practical: in the highest resolution mode (37.5 m 
range gate spacing), the A/D converters on each ACQ board acquire 12-bit samples sign- 
extended to 16 bits at a rate of 4 MHz. The conversion of three logarithmic and four 
linear channels (seven total) results a raw-data bandwidth of 56 Mbytes/second. At this 
rate, an eight-hour flight will necessitate storage of 1.6 Terabytes of information! By 
performing real-time reflectivity and Doppler processing and storing only the products, the 
data rate can be reduced by a factor of 400 or more (depending on the chosen integration 
cycles) allowing for convenient data archival on a 5 Gb Exabyte tape. A block diagram is 
provided in Figure 3-5 to illustrate the flow of information between the data system and 
other system components. The individual circuit boards that comprise the data system are 
described in limited detail in the remaining paragraphs of this section. 
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ER-2 DOPPLER RADAR DATA SYSTEM 













3.2.1 Radar Interface Board (RIB) 


The Radar Interface Board (RIB) is responsible for generating the various system clocks 
and phase-locking them to the 10 MHz STALO oscillator in the radar transmitter. It is 
further used to control radar transmitter and receiver operations by supplying a transmit 
waveform used to modulate the RF carrier, selecting receiver bandwidth, and performing 
radar calibration operations. Automatic gain control is accomplished by manipulating RIB 
supplied analog signals to 60 dB RF- and IF-attenuators in the receiver. Acquisition of 
analog system status information such as temperature and pressure measurements within 
the transmitter, receiver, and data system, is achieved using a 16-channel 12-bit A/D 
converter. Furthermore, an octal 8-bit DAC provides the board with analog output for 

signaling and control. This capability is enhanced by 16 bits of digital I/O which are not 
presently utilized. 

Central to the operation of the RIB is a single programmable AT&T WEDSP32C digital- 
signal-processor chip which executes embedded software to provide radar control. A 
typical embedded program appears in Appendix A and was used during the 1993 CAMEX 
experiments at the NASA Wallops Flight Facility. 

This program begins by assigning pointers to the Board Control Register (BCR) and 
Radar Status Register (RSR). These registers, shown in Tables 3-2 and 3-3, specify the 
radar operating parameters (e.g. PRF, pulse width, etc.) and indicate the operating state of 
the system. Immediately after system powerup, the radar enters a warmup mode, heating 
the TWT to the proper operating temperature before transmitting. During this warmup 
period the RSR reports an operational status of STANDBY, and the RIB embedded code 
acquires and records the following status information twice each minute: 
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• Board Status Register 

• Radar Status Register 

• Transmitter Base Temperature 

• Receiver Base Temperature 

• Transmitter Air Temperature 

• Receiver Air Temperature 

• Transmitter Local Oscillator Temperature (ADCDAC4) 

• Unused Channel (ADCDAC5) 

• Transmitter Nitrogen Pressure (ADCDAC6) 

• Unused Channel (ADCDAC7) 

• RIB Ambient Temperature (ADCDAC8) 


(BSR) 

(RSR) 

(ADCDACO) 

(ADCDAC1) 

(ADCDAC2) 

(ADCDAC3) 


The ADCDAC registers are 12-bit sampled analog signals connected to various 
thermistors and pressure transducers. 


After the pilot moves the operation switch to the RADIATE position and warmup is 
complete, the transmitter attempts to enter a radiate mode. To facilitate a successful 
mode change the RIB must immediately provide a trigger signal by asserting PRFSTART 
within the Board Control Register. If the RIB fails to provide this signal, and the 
transmitter is otherwise working properly, the transmitter will assert RSR[0] to indicate a 
first-level fault (FAULT 1). This kind of fault (first-level) will also be generated if sensors 
indicate the transmitter or receiver are not working properly. For example, the following 
conditions will all initiate a first-level fault: transmitter or receiver outside of temperature 
operating range, TWT high-voltage off, transmitter power low , or transmitter VSWR out 
of range. If no faults are indicated, the transmitter will enter its normal operational state. 
However, if some failure results in a FAULT 1 state, the transmitter will attempt to reset 
itself and enter the RADIATE state a maximum of three times. If the transition remains 
unsuccessful after the three attempts, the transmitter will completely shut down, assert 

RSR[1] to enter a catastrophic failure or second-level fault (FAULT2), and remain in this 
state until power is cycled. 
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Bit 

Net Name 

Destination 

Pin# 

Function 

0 

RSRDEBUGN 

PBCTL2 PAL 

11 

Radar Status Register debug 10) 

1 

PIOROUTN 

PBCTL2 PAL 

10 

Programmable IO are Outputs (0) & Inputs (1) 

2 

FREQINP 

PLL Clk Driver 

6 

Freq I/P Select b/w STALO (1) & Local Osc 10) 

3 

ADCDIFF 

INAIS 

1 

Select Differential f 1) or Non-diff 10) on ADCs 

4 

PRFCTL 

PRFCTL PAL 

El 

Select PRF in combination with PRF bit - 
HI: 8.8/1.1 KHz LO: 4.4/2.2 KHz 

5 

PRFSTART 

PRFCTL PAL 

E2 

Start PRF 

6 

OPTPULW 

PRFCTL PAL 

B6 

Select Pulse Width in combination with the 
PULSWID bit - HI: 0 . 54 s LO: 1.0/0.25 4 s 

7 

AGCSTART 

AGCCTL PAL 

2 

Start AGC 

8 

DSPINTO 

DSPINT PAL 

2 

DSP Interrupt Control Bit 0 

9 

DSPINT1 

DSPINT PAL 

3 

DSP Interrupt Control Bit 1 ~j 

10 

DSPINT2 

DSPINT PAL 

4 

DSP Interrupt Control Bit 2 

11 

DSPINT3 

DSPINT PAL 

5 

DSP Intemipt Control Bit 3 

12 

LDPCTL 

7-SEG DISPLAY 

4 

Left Decimal Point Control 

13 

RDPCTL 

7-SEG DISPLAY 

10 

Right Decimal Point Control 

14 

DONE 

VME CTL2 PAL 

11 

Signal End of Acquisition 

15 

NEWDWELL 

AGCCTL PAL 

7 

New Dwell Signal from DSP 


Table 3-3: Radar Status Register - RSRMS-n 
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During the normal operational state RADIATE, the RJB records status information 
frequently (-3 second intervals) and monitors RSR[9] for assurance that all is well. An 
RSR[9] status change from RADIATE to STANDBY might indicate a transmitter or 
receiver first-level fault due to some random glitch; however, it could also indicate the 
pilot has turned off the transmitter! To determine which is the case, the RIB waits ten 
seconds. If the event was a transient fault, the minor glitches will work themselves out 
during this time and the radar will automatically return to the RADIATE state. If the 
system remains in STANDBY for more than ten seconds, it is assumed that the pilot has 
turned off the transmitter or that a catastrophic second-level fault has occurred. 

Regardless, he RIB deasserts the tngger and signals the E-6 host computer to terminate 
acquisition by closing files. 

The RIB DSP32C processor receives a hardware interrupt each time the transmitter fires a 
pulse. These regular interrupts form the basis of internal program timing and provide 
synchronization among all processors. Every five minutes, the transmitter is switched into 
a calibration mode by asserting CALIB of Radar Control Register Zero (RCR0[3J). 
Attenuation settings are then read from a table and applied to the calibration attenuators 
by writing to the Radar Control Register One (RCR1). Presently, those calibration values 
step from 0 dB to -120 dB in 10 dB increments. 
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Bit 

Net Name 

Pest 

° 

Pin # 

Function 

0 

RSRVD 




1 

PULSWID 

J1 

11 

Pulse Width - HI: 0.25 us LO: 1. 0/0.5 us 

2 

PRF 

J1 

12 

Pulse Repetition Freq - HI: 2.2/1. 1 KHz LO: 4 4/8 8 KH 7 

3 

CALIB 

J1 

13 

Calibrate Mode - HI: Off LO: On 

4 

RECBW 

J1 

14 

Receiver Bandwidth - HI: 2.2 MHz LO: 8 MHz 

5 

RSRVD 




6 

PHASESHO 

J2 

1 

Phase Shifter Bit 0 - HI: 5.6° LO: 0° 

7 

PHASESH1 

J2 

— 

2 

Phase Shifter Bit 1 - HI: 1 1.3° LO: 0° 

8 

PHASES H2 

J2 

3 

Phase Shifter Bit 2 - HI: 22.5° LO: 0° 

9 

PHASESH3 

J2 

4 

Phase Shifter Bit 3 - HI: 45.0° LO: 0° 

10 

PHASESH4 

J2 

5 

Phase Shifter Bit 4 - HI: 90.0° LO: 0° 

11 

PHASESH5 

J2 

6 

Phase Shifter Bit 5 - HI: 180° LO: 0 ° 

12 

PHASESHLAT 

J2 

7 

Phase Shifter Latch - HI: 0=D LO: 0=0 

13 

RSRVD 



— ^ 

14 

RSRVD I 




15 

RSRVD 





Bit 

Net Name 

Dest 

Pin# 

Function 

0 

ATNAO 

J2 

8 

Attenuator A Bit 0 - HI: 1 dB LO: 0 dB 

1 

ATNA1 

J2 

9 

Attenuator A Bit 1 - HI: 2 dB LO: 0 dB 

2 

ATNA2 

J2 

10 

Attenuator A Bit 2 - HI: 4 dB LO: OdB f 

SB 

ATNA3 

J2 

11 

Attenuator A Bit 3 - HI: 8 dB LO: OdB 

IS 

ATNA4 

J2 

12 

Attenuator A Bit 4 - HI: 16dB LO: OdB 

5 

ATNA5 

J2 

13 

Attenuator A Bit 5 - HI: 32 dB LO: 0 dB 

6 


J2 

14 

Attenuator B Bit 0 - HI: 1 dB LO: 0 dB 

7 

ATNB1 

J2 

15 



Attenuator B Bit 1 - HI: 2 dB LO: 0 dB 

8 

IATNB2 

J2 

16 

Attenuator B Bit 2 - HI: 4dB LO: OdB 

9 

ATNB3 

J2 

17 

Attenuator B Bit 3 - HI: 8 dB LO: OdB 

10 

ATNB4 

J2 

18 

Attenuator B Bit 4 - HI: 16dB LO:OdB 

ii 

ATNB5 

J2 

19 

Attenuator B Bit 5 - HI: 32 dB LO: 0 dB 

12 

RSRVD 




13 

RSRVD 




14 

RSRVD 




15 

RSRVD 
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Table 3-6: Board Status Register - BSR[15:01 


Net Name 


VMEQBAFN 

VMEQBAEN 

DATDQAAFN 

PATDQAAEN 

DATDQBAFN 

PATDQBAEN 

DSPINQAFFN 

DSPENQAEFN 


Source 


VMEQB 

DATDQA 

PATDQA 

DATDQB 

DATDQB 

DSPINQA 

DSPINQA 


Functi on 

VME 

VME Queue B Almost Empty Fla f 

Data Destination Queue A Almost Full Fla p 
Data Desttination Queue A Almost Empty F lat 
Data Desttination Queue B Almost Full Fla f 
Data Desttination Queue B Almos t Empty Fla* 
DSPINQA Full Flag 
DSPINQA Empty Flap 


jFFN 
NGAINQEFN 
QGAINQFFFN 
QGAINQ EFN 
RNRNA 
RNRNB 


3.2.2 Data Acquisition Board (ACQ) 

Received echoes are sampled by dual 4 MHz four-channel 9U VME boards featuring 12 

bits of resolution per sample. Each board operates under the direction of a DSP32C 

digital signal processor chip whose job is to compute destination tags for each sample 

which are then paired with the samples acquired by the A/D's as routing information. The 

High-Speed Data Bus Controller (HSDBC) uses these tags to route each sample to some 

desired location, often a processing-zone on the LMAP Unit Array Processor Board 

(LUA-200), much like the Post Office uses address information on a letter to route it to 
the desired house. 


The software used to scatter EDOP raw data for reflectivity processing is provided in 
Listing 2. A hardware interrupt notifies the program each time a pulse has been 
transmitted. After each interrupt, the software queues up destination tags for all 436 
range gates to be acquired this pulse. This is accomplished by reading a table of 
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predetermined addresses. Four processors are required for reflectivity processing; hence, 
the first 109 tags route samples to Zone 0, the next 109 are routed to Zone 1, and so on. 
After sufficient data for an entire dwell have been scattered for processing, the global 

token (see §3.2.3) is passed to the LUA so that its results can be delivered from its output 
queue. 

3.2.3 LMAP Unit Array Processor Board (LUA-200) 

The LUA-200 is a highly-parallel signal-processing board containing a Global I/O 
Interface, Local I/O Interface, VME Interface, Input and Output Data Queues, and eight 
linearly-mapped computational zones. The reader is referred to Figure 3-6 for a block 
diagram explaining the relationship between these components. 

Raw data arrive at the board for processing via two 40-bit bi-directional high-speed data 
buses (HSDB) dedicated to the transfer of information between the ACQ and LUA 
boards. Each bus is comprised of 32 data lines, 4 control lines, 4 lines for destination tags, 
and runs at !4 the system clock (10 MHz) to achieve a maximum throughput of 50 
Mbyte/sec/bus. 


Inter-zone communications is facilitated by hardware-based Input and Output data queues 
Each queue is a 40-bit x 5 12- word FIFO possessing a global and local side. The global 
side of the queue is connected to the Global I/O Interface which manages the transfer of 
information between the I/O queues and the HSDB. The Local I/O interface similarly 
manages information transfer between the I/O queues and processing zones. Data (i.e. a 
destination tag/data word pair) received from the HSDB by the Global Interface controller 
are placed in the Input Data Queue of the appropriate bus where they remain until their 
destined zone requests that they be read from the queue and transferred to the zone's local 
memory. The data are removed from the queue and flexibly routed to the desired zone 


36 




MARTIN MARIETTA 
AERO A NAVAL SYSTEMS 
LUA200 : IMAP TOP LEVEL~ BLOCK DIAGRAM 
D I 1 M K. PLANTE 
DATE : 2/1/91 I SHEET 01 CF XX 





















through a series of cross-point data switches on each board — the Local I/O Interface. 
Because the queue is based on a FIFO architecture, they are read out in the same order 
they are received. If a processing zone requests additional data and data destined for that 
zone are not at the front of the queue, the zone will block or wait until such data are 
available. We say that the input data stream is flow controlled. 

When a zone is finished processing and desires to transfer information to another location, 
it does so by writing the data and appropriate destination tag to the Output Data Queue. 
To control the order is which output data are queued up, only one zone is given 
permission to write at any given time. This permission or write arbitration is governed by 
the use of a semaphore called the local token , and there is one local token for each local 
bus. Any zone attempting to access the Output Data queue without possessing the write 
access privileges granted by the local token will be held in a pending-access state until the 
zone receives the local token from the current token holder. 


After data are placed in the Output Data Queue, the Global I/O Interface reads the output 
destination queue and passes its associated data word to the proper destination. If the 
data are destined to another zone on the same board, the Global I/O Circuit merely places 
them back in the Input Data Queue. If they are instead destined for another board, the 
Global I/O Circuit checks to see if the board is a global token holder. The global token 
provides bus arbitration on the HSDB in a manner similar to the local bus so that only one 
board in the system has access to the HSDB at any given time. If the board is possesses 
the global token, the Global I/O interface removes data from the Output queue and passes 
them to the HSDB Controller for transfer to the Global I/O Controller of the destination 
board (via the high-speed data bus). If the board does not possess the global token, the 

data remain in the Output queue until the global token is received by the current token 
holder. 


38 


Three operating modes for the LUA Board, the Master, Slave, and Local Access modes, 
further govern Global I/O. 


Master Mode allows write access to the High-speed Data Bus and is allowed only by the 
Global Token holder. In this mode, the presence of data in the Output Data Queue 
initiates a data transfer to the destination specified in the Data Destination Queue. Valid 
destinations include the Data Input Queue of the same board. Data Input Queue of a 
different board, or any global bus command (Figure 3-7). 

Boards which do not possess the Global Token are considered to be in Slave mode and 
will remain Slaves until the Global Token is relinquished by the Master. A Slave can 
receive data from its Data Input Queue and additionally write to the Data Output Queue 
as Local Token issues permit. However, no data will be read from the Output Queue until 
the board becomes a Master by receiving the Global Token from the current Master. 

In some cases it is desired to pass data only between zones of the same board. Because no 
data are transferred over the HSDB, it does not make sense to restrict data flow with a 
global token; however, a the Global I/O controller will not read the Output queue of a 
Slave. To address this problem, a third mode. Local Access mode, is available. In Local 
Access Mode, the board isolates itself from the HSDB allowing data to be passed through 
the Output Queue back to the Input Queue without ever possessing the Global Token. In 
this mode, the Global DMA Controller processes data from the Data Output Queue as 
though the board was a Master— except that destinations to other boards are invalid. Any 
data destined to different board will be lost. Similarly, any data transferred to a board in 
Local Access mode over the HSDB will be lost. To enter Local Access Mode, the board 
must first be the Global Token holder (i.e. Master). After setting HCR[1], the board must 
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DESTINATION 

BOARD 

CODE 

NODE 

(HEX) 

(DEC) 


DSP 

ZONE 

(DEC) 


REGISTER! 

FUNCTION 


0(000) 

7(007) 



DATA 

DATA 


MASK REG 


CTLR REG 


NODE 


DATA 


MASK REG 


TOKEN 


DATA 


MASK REG 


TOKEN 


DATA 

DATA 


MASK REG 


CTLR REG 


NODE 


DATA 


MASK REG 


TOKEN 


ZONE VME 

OPERATION OPERATION 


WR/RD 

WR/RD 


WR/RD 


WR/RD 


WR/RD 



BCAST WR 


BCAST WR 


PASS 


BCAST WR 


BCAST WR 


PASS 


BCAST WR 


BCAST WR 


PASS 


BCAST WR 


BCAST WR 


PASS 


DATA 

WR 

WR/RD 

DATA 

WR 

WR/RD 

MASK REG 

WR 

WR/RD 

CTLR REG 

WR 

WR/RD 

NODE 

WR 

WR/RD 



BCAST WR 


BCAST WR 


PASS 


WR/RD 

WR/RD 


WR/RD 


WR/RD 


WR/RD 


BCAST WR 


BCAST WR 


PASS 


Figure 3 7: Data Destination Codes (courtesy of Martin Marietta) 










































































































pass away the Global Token. The board will remain irj Local Access Mode until resetting 
HCR[1] at which time it will become a Slave. 

Each of the eight linearly-mapped parallel processing zones is centered around an AT&T 
WEDSP32C digital signal processor housed in a standard 133-pin square pin-grid array 
(PGA) package. CMOS fabrication technology provides high-speed operation with low 
power consumption and two execution units increase computational performance. The 
Data Arithmetic Unit (DAU) utilizes four 40-bit floating-point accumulators to store 
operation results and data type conversions. A 32-bit floating-point adder and 40-bit 
floating-point multiplier work in tandem to jointly perform 25 million floating-point 
multiply-accumulate operations per second (peak) when clocked at the maximum speed of 
50 MHz. Integer operations, as well as control, logical, and data move operations are 
conducted by a second execution unit — the Control Arithmetic Unit (CAU) — at a rate of 
12.5 million instructions per second (MIPS). Memory internal to the processor 
includes 1536 32-bit words allocated in three 512-word banks. An external memory of 
32K x 32-bit 20-ns zero wait-state static RAM per zone provides ample room for program 
and data storage. Three 40-MHz LUA-200's are utilized by the EDOP data system 
yielding an aggregate 480 MFLOPS of computational power for Doppler extraction, 
reflectivity averaging, and automatic gain control. 

3.2.4 Host Computer (E-6) 

An E-6 Eltec single-board computer based on the 68030 microprocessor is utilized as a 
host. This host is responsible for initializing various system registers at power up and 
loading all DSP32Cs with software stored on a battery-backed non-volatile 8 Mb RAM 
drive. Status information are additionally stored on this RAM drive by the E-6 which is 
also responsible for archiving processed data sent to it via the VME bus during operation. 
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A 1.2 Gb Micropolis hard drive is currently used to stpre processed results, and a 5 Gb 
8505C Exabyte streaming tape is being evaluated for future use. 
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Chapter 4 - Algorithm Development and Analysis 

4.1 Reflectivity Processing 

Conventional weather radars relate the received echo power from a precipitating target to 
its radar reflectivity factor, Z, through a special form of the radar equation applicable to 
distributed targets. The radar reflectivity factor can then be empirically related to the rain 
rate. Nevertheless, the problem of accurately assessing Z based on power measurements 
is not without difficulties. As wind turbulence shuffles the hydrometeors, backscattered 
microwaves undergo constructive and destructive interference causing the magnitude of 
received power to vary greatly. To effectively estimate the true mean received power, one 
must average many power samples to reduce the variance of this fading process. 

4.1.1 The Need for Reflectivity Averaging 

To better illustrate the nature of fading, let the receiver input voltage due to the i* 
scatterer be expressed in polar form as: 

V ie * 

where V; is its magnitude and tfc is the instantaneous phase. We are normally interested in 
the amplitude and phase of the complex voltage, V, due to a collection of scatterers and 
can easily express this as a superposition of received voltages from individual scatterers: 

V = IV,e*' 

i=i 
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Of course, the sum of a complex-number series is itself a complex number, so it is 
convenient to express the complex voltage, V, in terms of an envelope, V e , and a phase 
angle <)>: 

V = V e JI> 


If the number of randomly phased scatterers is sufficiently large, the central limit theorem 
provides reasonable assumption that the received voltage is a bivariate Gaussian 
distribution of the form 


p(V e ,<|>) 


Ve 

2n a 2 


e -vi/ 2 CT ; 


The individual density functions for V e and <h can then be determined by "integrating out" 
the opposing variable. 

p(V e ) = Tp(Ve.4>)d4> = — e - v ' /2 ° ! 
o <r 


pf^) = Jp(V e ,4 ) )dV e =—• 
o 271 


It is interesting to observe that p(V e ) is the familiar Rayleigh distribution whose first two 
moments are given by 



V? = 2a 2 
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The variance can then be determined as the difference between the second moment and the 
square of the first: 

<7v t = Y? - ' Ve = (2 - J j a 2 = 0. 429 a 2 

Assume the power in the envelope voltage is developed across a 1-ohm resistor. The 
envelope power can be expressed as the square of the envelope voltage. 

P=V 2 

If we further note that the differential power is related to the differential voltage as 
dP = 2 V e d V e 

then the Rayleigh distributed envelope voltage can be converted to a distribution for 

power by transforming the voltage variable to obtain an expression for the distribution of 
power at the radar receiver input 

P(P)dP = p(V e )dV = — e' p/2 ° 2 


Notice that the power is distributed exponentially whereas the voltage was Rayleigh 

distributed! The first two moments of the exponential power distribution are easily seen 
to be: 


P = 2 a 2 

P 1 =2p 1 
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revealing a classic property of exponential distribution? wherein the standard deviation is 
equal to the mean! 

ap = P 2 - p =p* 

(Tp= P 

By now it should be evident that the power fluctuations due to fading can be significant 
and some type of estimator based on multiple observations must be employed to reduce 
this variance to acceptable levels. Before we can address that issue, however, we must 
follow the received pulse through the receiver to understand its effects on the fading 
statistics. Weather targets, especially convective targets, exhibit a range of echo strengths 
exceeding 80 dB. This rather large variation places tough demands on the receiver design 
and frequently logarithmic detectors, whose output signal is proportional to the logarithm 
of the input power, are used to minimize receiver saturation. Let L be defined as input to 
the log-receiver, expressed in logarithmic form, where P is the instantaneous, 
exponentially-distributed received power: 

L = 1 0 log(P/p, ) 

The reference value P, is chosen to be 1 mW so that the units of power-level L are dBm. 
The mean power can be similarly written: 

Lo = 101og(P/p,) 

Certainly, the logarithmic transfer function compresses the range of output, but 
additionally has less intuitive effect of adding a statistical bias to the measurement. To see 
that this is true, it is again necessary to transform the power variable using 
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L-L o = 101og(P/P) 


As might be expected, the distribution of log-power is no longer exponential but is instead 
a slight variation given by (Sauvageot, 19xx) 

P(L) = m exp[m(L - L 0 ) - e m(L_u) ] 

with m = ln(10)/10. The first two moments, can be obtained by integrating the density 
function 

E{L}= jLmexp[rn(L-Lo)-e m(L ' u>, ]dL 
E{L 2 } = j L 2 m exp[m(L - L 0 ) - e m(L_u) ]dL 


These equations were numerically integrated to determine the mean and variance of the 
detected spectrum. By definition, the true power mean is L 0 . However, the expected 
value of L is given by E { L } = L 0 -2. 51, indicating that the estimate is statistically biased. 

This means that one must add 2.51 dB to the estimate to get the true mean. The variance 
is nearly constant at 3 1 dB. 

4.1.2 An Averaging Algorithm for the DSP32C 

Clearly, the power estimate must be improved if it is to be scientifically useful. A 
particular type of estimator known as the maximum-likelihood estimator was chosen to 
provide mean-power estimates for the EDOP application. It works as follows: Let 

/(x,,...,x n l0) be the joint probability density function of the random variables X, ... X. 

where 6 is the parameter to be estimated— in this case, the mean power. This function is 
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called the likelihood function and represents the likelihood or certainty that the values x,, 
x 2 , .... x n will be observed when 0 is the true value of the unknown power mean. It is then 
reasonable to assume that 0 can be estimated by choosing its value such as to maximize 
the likelihood of observing x„ x 2 , x n . In other words, the maximum-likelihood estimate 
is defined to be that value of 0 which maximizes the likelihood function. 

The maximum-likelihood estimator for the mean of an exponential distribution is given by: 


To generate estimates of the mean echo-power received by EDOP, a real-time embedded, 
parallel algorithm was developed for the EDOP data system to average a predefined 

number of samples, N, and archive the results. A listing of this software is included in 
Appendix A. 

To implement the algorithm, two constraints were observed to maximize the performance: 

1 . Keep the pipeline full 

2. Minimize use of memory 

As previously mentioned, the DSP32C is a pipelined digital signal processor. Pipelined 
architectures execute their instructions in assembly-line fashion to realize superior 
performance. Each stage of the pipeline completes a small part of the instruction being 
executed in a fraction of the time needed to complete execution of the entire instruction. 
Moreover, each stage is connected together to form a pipe so that instructions enter the 
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pipe at one end, undergo various stages of execution in each pipe segment, and exit other 
end. 

Throughput of the pipeline is the rate at which completely executed instructions exit the 
pipeline. Because multiple instructions are overlapped in execution, an instruction can exit 
the pipe with each machine cycle resulting an increase in CPU performance direcdy 
proportional to the length of the pipe. However, superior performance is sustained only 
so long as the pipeline is kept full. If the pipeline is not kept full, some segments will be 
idle when they could be executing instructions and the throughput will decrease. As this 
bubble of inactivity propagates through the pipe and exits the other end, one or more 

machine cycles will be wasted while the CPU waits for the next instruction to exit the 
pipe. 

The data arithmetic unit (DAU) of the DSP32C employs a four-stage pipeline to perform 
25 million floating-point operations per second. Its instruction set revolves around a 
multiply-accumulate operation wherein a floating-point multiplier and adder work in 
parallel to perform computations of the form a = b + c*d. Each DAU multiply-accumulate 
operation involves three floating-point operands. Two of these are multiplied together 
and added to the third. The final result is stored in an accumulator as an intermediate 

result to be used in subsequent operations and can be additionally stored in memory or 
sent to an I/O unit. 

The DAU supports several multiply-accumulate instruction formats. However, for the 
purpose of discussing the DSP32C pipeline, only one format will be examined: 

Z = aN = aM + Y*X 
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This instruction is executed in four stages as follows: 

1. XY fetch 

2. multiply Y * X 

3. accumulate product with aM 

4. write result to memory (optional) 


If several multiply-accumulate operations are executed back to back, the DSP32C fills the 
pipeline such that one stage of each instruction is completed every machine cycle. For 
example, consider the following block of multiply-accumulate instructions: 

1. aN = aM + Y, + X, 

2. aN = aM + Y 2 + X 2 

3. aN = aM + Y 3 + X 3 

4. aN = aM + Y 4 + X 4 

5. aN = aM + Y 5 + X 5 


where X, X 5 and Y,— Y 5 have been subscripted to aid visualization of the data flow 
through the pipeline. During each machine cycle, the four instructions will be in various 
stages of simultaneous execution within the pipeline: 


Table 4-1: Pipeline Data Flow 


Cycle 

Segment 4 

Segment 3 

Segment 2 

Segment 1 

1 




XY fetch, 

2 



multiply, 

XY fetch. 

3 


accumulate, 

multiply. 

XY fetch. 

4 

write, 

accumulate. 

multiply. 

XY fetch,, 

5 

write. 

accumulate^ 

multiply,, 

XY fetch. 

6 

write. 

accumulate. 

multiply. 

.,.1 

7 

write 4 

accumulate. 



8 

write. 
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It should be evident from the table that the intermediate accumulated result from the first 
instruction will not be available for use in subsequent multiplier operations until Cycle 4. 
Hence there is a latency of three instructions before the intermediate results from the first 
instruction can be used. In other words, if it is desired that the intermediate results of one 
instruction be used in a subsequent instruction, no less than three unrelated DAU 
instructions must be executed before the intermediate results are available for use. Such 
latencies are inherent in pipelined architectures and the DSP programmer must design 
software carefully, keeping the pipeline full to maximize the throughput. 

To begin the actual coding of software for this digital-signal processor, a slightly different 
approach is used than for general purpose CPUs. In this approach, the DAU or floating- 
point unit is considered to be the central processing unit of the DSP32C while the CAU is 
viewed as a co-processor to generate addresses and perform occasional logical operations. 
To fully utilize the processing power of the parallel array of DSP32Cs, consideration must 
be given to how the processing will be split among the processors while maximizing the 
pipeline efficiency. For the reflectivity averaging software, it was decided that each 

processor would perform the same estimation algorithm on a subset of the range gates to 
be processed. 

After each RF pulse is transmitted, three A/D's provide range-gated samples of received 
power on each of three channels — nadir antenna, forward co-polarized antenna, and 
forward cross-polarized antenna . These echo power measurements are sent to the data 
system via both 32-bit HSDB buses, with two 16-bit channels packed into a single 32-bit 
word per bus. The reflectivity program expects that each word read from HSDB FIFO A 
contains channel one in the most-significant 16 bits while channel two is stored in the 
least-significant 16 bits. Similarly, channel three is expected to reside in the most- 
significant 16 bits of FIFO B; however, channel four is ignored. A typical flight 
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configuration is to place the Nadir Co-pol. on channel.one, the Forward Co-pol on 
channel two, and the Forward Cross-pol on channel 3. These range-gated measurements 
are placed into a buffer for accumulation by the following zero-overhead core loop of 
reflectivity DSP32C-assembler code: 

do 5, rl 8 

al = float(*rlO++) 
aO = float(*rlO++) 
a2 = float(*rll+rl5) 

*r6++ = a3 = *r6 + al 
*r5++ = a3 = *r5 + aO 
*r7++ = a3 = *r7 + a2 

The first line indicates that the following six (5+1) lines of code are to be executed N 
times, where N is the value contained within hardware register rl8 at the time of 
execution. In the buffer, addressed by hardware pointers rlO and rl 1, the received-power 
samples remain packed together with two channels per word. The next three instructions 
in the loop extract the 16-bit integer samples from the buffer and convert them to a 32-bit 
DSP32C proprietary floating-point format, placing each result in a 40-bit accumulator. In 
this example, samples from channel one are placed in accumulator al while samples from 
channels two and three are placed in accumulators aO and a2, respectively. The remaining 
three instructions accumulate the measurements into separate range-gated tables addressed 
by hardware pointers r5, r6, and r7. At the end of each integration period, the 
accumulated results are placed in the Output queue along with appropriate destination 
tags sending the results to the VME bus for subsequent archival by the E-6 host. A pass 
global-token command is then queued to return the global token to the ACQ board so the 
HSDB Controller can deliver additional samples to the LUA for processing. 
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4.1.3 Number of Independent Samples 

The reflectivity algorithm used by EDOP averages a series of observations, Lj, for one 
second to reduce the variance to acceptable levels. However, the observations are 
partially correlated, so to determine how much the variance is reduced through averaging, 
one cannot simply reduce the variance by the number of averaged samples order. Instead, 
the variance is reduced by the number of effectively independent samples used in forming 
the average. Traditionally, samples are considered to be independent when the correlation 
between samples approaches zero. 

Several authors (e.g. Marshall and Hitschfeld, 1953; Sauvageot; Walker et. al, 1980) have 
attempted to develop receiver models for the purpose of establishing the rate of sample 
decorrelation. For example, Kerr (1951) provides the following relation between the 
normalized input and output correlation functions of a logarithmic receiver: 

P^W = 4^P 1 ”m' ! 

K m = l 


Here ’ Piog^) is the desired normalized output correlation and p(x) is the correlation 
function of the input process. If the input signal is assumed to be Gaussian correlated, 

p(x) = exp(-t 2 /2 a 2 ), 

then the output signal is correlated as 

Piog(^) = ~r Xm' 2 exp(-mx 2 /a?) 

It m = l 
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One definition for the number of independent samples js given as a ratio of the variance of 
a single sample to the variance of the sampled mean (Lee, 1961; Nathanson, 1969). For 
discrete sampling, (Walker et. al, 1980) give 


Nr ' - 77 + 77 CN q ~ K q ) p q (K q l q ) 


2 N,-l 


N, 


Nq K,=l 


where N q is the number of samples spaced l q apart where q represents time. Thus to 

determine the number of independent samples, N it one merely substitutes the appropriate 
correlation function. 

These efforts are of somewhat questionable utility because it is so difficult to specify the 
meteorological processes that are viewed by the airborne receiver. Highly turbulent 
conditions will certainly cause rapid decorrelation between samples and the ensuing large 
number of independent samples will result excellent mean-power estimates for even short 
integration periods. On the other hand, stratiform conditions with weak air motions might 
push decorrelation times on the order of seconds so that observed sample decorrelation 
results primarily from aircraft displacement. In view of the consequence of incorrectly 
specifying the observable meteorology, it is best to be conservative. 

For airborne radars, it is generally accepted that pulses are decorrelated in the time it takes 
the aircraft to move one beamwidth. The EDOP radar flies with 0.76 m antennas on an 
ER-2 platform possessing an airspeed of approximately 200 m/s; therefore, pulses are 
expected to decorrelate in ~4 ms due to aircraft motion, providing a minimum of 250 
independent samples during a 1 -second integration period. Again, certain meteorological 
conditions (e.g. severe convective regions of a storm) can decorrelate the pulses at a faster 
rate. It seems the best way to determine independent samples obtained on a particular 
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flight is to examine the taken measurements. If the conservative decorrelation estimate of 
4 ms is used, the expected post-integration ( 1 second) variance is 

_2 

a ‘ vg = 250 =31dB ~ 101 ° g(250) “ 702dB 

4.2 Doppler Processing 

The EDOP radar is primarily intended to study deep, well-developed, precipitating 
systems where both horizontal- and vertical-wind velocities can exceed 40 m/s. The 
desired accuracy of velocity estimates is on the order of 0.5- 1.0 m/s. However, errors 
arising from signal-processing operations and aircraft platform instabilities can induce 
significant errors in the vertical velocity that are sometimes comparable to the 
meteorology of interest (Heymsfield, 1989). The principle sources of velocity 
measurement error for the EDOP system result from uncertainties in the auto-covariance 
estimator and antenna pointing errors as discussed in §4.2.3 and §5.2. 

4.2.1 Autocovariance Processing 

The autocovariance C xx (t,,t 2 ) partially describes the time-domain structure of a random 
process X(t) and is defined: 

Cxx(ti,t2) = Rxx (t! , t 2 ) — M^x (ti)M- x (t2> 

where the autocorrelation is given by the expected value of the product of the process 
at two times, t, and t,: 

Rxx(t,.t 2 )-E{X'(t,)X(t 2 )} 
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Loosely speaking, a process is called stationary if its distribution functions or certain 
expected values are invariant with a shift of the time axis [Shanmugan, 1988], If the 
autocorrelation function depends only on the time difference t but is not explicitly a 

function of times tj and t,, and the mean is constant, then the process is said to be wide- 
sense stationary. 

R xx( T )= E{X*(t)X(t + t)}, 

E{X(t)} = m = constant, 

The autocorrelation function of a stationary random process tells us something about how 
rapidly we can expect the random signal X(t) to vary as a function of time. If the 
autocorrelation function rapidly decays to zero, we can expect the signal X(t) to vary 
rapidly with time. If the autocorrelation decays slowly to zero, then X(t) will be a slowly 
changing process. Likewise, a process X(t) with periodic components will exhibit a 
periodic autocorrelation. Thus, we may correctly conclude that the autocorrelation 
function provides some information describing the underlying spectral content of the 
random process. The exact relation is given by the Wiener-Khinchine relationship which 
states that the frequency-domain power distribution or power-spectral density of a random 
process X(t) is given by the Fourier transform of its autocorrelation function: 


Sxx(f) - F{R xx (t)} = jR xx (T)e' j2 " f, dt 


For any given range gate, the receiver produces a series of noise-corrupted complex 
voltage samples, V(kT s ), separated by the interpulse period, T s Each sample can be 
written as the sum of signal and noise: 

v (kT.) = Vk + nk, k = 0,l,...,M-l 
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The signal term can be alternatively expressed as a Doppler shifted weather-echo signal 
whose spectrum is centered about zero frequency: 

V (k T s ) = Sk e ja)dkT ‘ + nk 


The autocorrelation function of the process is therefore 
R vv (mT s ) = E{V*(kT s )X[(k + m)T s ]} = Sp(mTs)e^ ,) ' ,mT, + N p 5 m 

where N p is defined to be the mean noise power, p is the normalized correlation function, 
and 5 m is 1 for m = 0 and zero elsewhere. The difficulty in evaluating this autocorrelation 
function centers around determining the normalized correlation function. An easier 
alternative is to determine the autocorrelation directly from the power spectral density by 
inverting the Wiener-Khinchine equation above: 

RxxOO = F-' {Sxx(f)} = jSxx(f)e j2,rf 'df 


For the purpose of determining statistical estimators, it is both convenient and reasonable 
to assume that the power spectrum has a Gaussian shape: 


S(v) = 



(u— l))M | 2 NpTs 

2 cl ) X 


Substitution yields 



Rw(mTs) = Sexp 


-8 


7CCT v mT! 


) je j4 ^ mT ^ + Np6 m . 


This autocorrelation function can then be compared with the one above to identify the 
normalized correlation function as: 


P(mT t ) = exp 


gf 7tq v mT s 

' l X 




At this point, it should start to become evident that if the complex autocorrelation function 
of the random process is known, the mean velocity can be obtained from its argument! 

arg[Rw(mTs)] = 47a)mT > , ms t0 

A. 


The simplest case is to solve for the mean velocity using a single time lag (m = 1 ). The 
resulting velocity estimate is then given by: 

\) = a/47iT s )arg[R vv (T s )] 


This well-known autocovariance or pulse-pair estimator is applied to each range-gate to 
provide estimates of the Doppler velocity. This algorithm was first described by Rummler 
(1968) and has since become the primary Doppler velocity estimator used by the 
meteorological community. A spectral- width estimator can also be derived as a function 
of the complex autocorrelation function and is based on the ratio of lag- 1 and lag-2 
estimates (Srivastava et al„ 1979): 


a 2 = 


247t 2 T, 2 


lnl 


R(T.) 


R(2T.) 
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The new variance estimator is of comparable quality to its more popular lag-l/lag-0 
counterpart (Benham et al„ 1972; Sirmans and Bumgarner, 1975) but does not explicitly 
require an estimate of the signal-to-noise ratio (SNR): 


(Told ' 


8tc 2 T? 


1 - 


jR(Ts)| 

R(0) 


(l + SNR-') 


Interestingly, the new variance estimator of Srivastava can be combined with the old one 
to provide an estimate of the SNR itself based on all three lags. The SNR can then be 
used as a statement of confidence for the reliability of these estimates: 

SNR = j: \MLil 

[R(0)|R(2T.) l/, |-|R(T.)r] 


4.2.2 DSP32C Implementation 

Doppler processing of the four linear channels represents the bulk of computations to be 
performed by the EDOP data system and centers around the computation of the complex 
autocorrelation function as described in §4.2. 1 : 

V ' = l2 i = 

Zj is the received complex echo-voltage of the i* pulse, X is the transmitted wavelength, N 
is the number of pulses per dwell, and T s is the interpulse period. 
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The autocorrelation is essentially a series of complex multiply-accumulate operations and 

is ideally handled by the DSP32C floating-point DSP. Consider the lag-1 autocorrelation 
for a dwell of N pulses: 


R(T S )~ IZpZ p+1 , 

P-0 


p = 0,l,...,N-l 


Both sides of the equation can be expressed as real and imaginary components: 


I(T,) + jQ(T,) = |[l p+J Q p ]-[ U + jQw ] 

= J; I 1 ' " ■* Q»I- *p +| + jQp+iJ 

N-l 

~ p ? 0 Ip Ip+I " J Ip+| Qp + Jip Q P+ i + Q P Q p+ i 


Collecting terms. 


I(TJ=I I p Ip+i+Q p Q p+ , 

p=0 

Q(T.)= I l p Q P+ i -l P+ i Q p 

p=0 H 


The lag-0 and lag-2 autocorrelation products are computed in a similar manner and 
summarized below: 


are 
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1(0) = i i;+q; 

p=o 

Q(o) = o 

I (Ts)= 2 Iplp +1 + Q P Q P+I 

p=0 

Q(T,)=S I t Q,„-l,. l Q t 

p=0 K 

I(2T s )= X IpI P+ 2 + QpQ p+2 

p = o 

Q(2T,)= Z I,Q«-IwQ p 

p~o 


To implement the algorithm in DSP32C assembler, it is again necessary to design the code 
so that the pipelined architecture of the CPU is optimally utilized. Shown below is the 
fourteen line core of EDOP's pulse-pair code (the complete listing is contained within 
Appendix B) which is used to accumulate autocorrelation products as described above. In 
flight, pulses will be correlated for a relatively long period of time, say, one-half second. 

At EDOP's maximum PRF of 4400 Hz, 2200 pulses will be correlated for each range gate. 
This loop correlates the 2200 pulses in small blocks to provide optimum utilization of the 

pipeline, accumulating the I and Q intermodulation products needed for the first three 
autocorrelation lags: 

do 12,rl6 
aO = *r6++ 
al = *r9++ 
nop 

a2 = *rl + aO * aO 
*rl = a2 = a2 + al * al 
a2 = *r2 + aO * *r7 
*r2 = a2 = a2 + al * *rl0 
a2 = *r4 - aO * *rl0++ 

*r4 = a2 = a2 + al * *r7++ 
a2 = *r3 + aO * *r8 
*r3 = a2 = a2 + al * *rll 
a2 = *r5-a0* *rll++ 


/* IpPrev */ 

/* QpPrev */ 

/* latency... */ 

/* a2 = 10 + Ip A 2 */ 

/* 10 = Ip A 2 + Qp A 2*1 
/* a2 = II + Iplp+l */ 
/*I1 =a2 + QpQp+l */ 
/* a2 = Q1 - IpQp+1 */ 
l* Q1 = a2 + Qplp+l */ 
/* a2 = 12 + IpIp+2 *1 
/* 12 = a2 + QpQp+2 */ 
/* a2 = Q2 - IpQp+2 */ 
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*r5 = a2 = a2 + al * *r8++ /* Q2 = a2 + QpIp+2 */ 

Hardware pointers r 1 — r5 point to range-gated tables where the accumulated 
intermodulation products are kept. Often the total number of pulses being processed is 
too large to simultaneously store them in memory. Therefore, a subset of these pulses is 
pulled into memory, processed, and then discarded so that an additional subset may be 
accessed. To ensure continuity of processing between the subsets, the last three pulses 
from each subset must be retained in memory as each new subset is created. These 
retained pulses are addressed by pointers r6 — rl 1. 

To determine the load (i.e. number of range gates and subset size) placed on each 
processor, simulations of the completed algorithm were run to calculate the number of 
instruction cycles needed to execute the algorithm for various load configurations. 

4.2.3 Pulse-pair Velocity Estimate Uncertainties 

An expression for the variance of the pulse-pair velocity estimate is given by Zmic (1977) 
for correlated but spaced pairs 

var('U) = A, 2 [327t 2 T 2 p 2 (T s )] * {M' 2 [l -p 2 (T,)] X 

M-1 

x xp (mT)(M-lml) + N 2 /MS 2 + (2N/MS)[1 + p(2T s )(l/M — 1 )§ t _ t 0 ]} 

where M is the number of sample pairs, and T is the spacing between pairs. If a Gaussian 
spectrum is assumed, this expression can be approximated by 

var(u) = X 2 [32tt 2 r 5 p 2 (T s )]'’ {[ 1 - p 2 (T,)] T s /2 a un T + N 2 /S 2 + 2(N/S) [1 + p(2 T,) 5t-t..o]} 
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For large signal-to-noise ratios, narrow spectrum widths, and contiguous pairs, this 
expression further reduces to 


var(b) « V(8MT s Vi) 


To illustrate the quality of this estimator, consider the following example based on a 3 cm 
radar. In a convective storm, a typical spectrum width of 5 m/s would result in a velocity 
estimate uncertainty of 0.047 m/s when integrating 1000 pulses at a PRF of 4400 Hz. 

These equations are derived using perturbation analysis. For very narrow spectrum widths 
or low signal-to-noise ratios, a large number of pulses must be integrated to obtain valid 
results. Specifically, the following two conditions are necessary to ensure validity of the 
estimated variances: 

2nM a„„ » 1, 

P 2 (T s )M » (N/S-M ) 2 

The first condition expresses the requirement for a large number of samples whereas the 
second condition ensures that the arg[R(T s )] is small compared with 2n. 
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Chapter 5 - Conclusions 


5.0 Status of Software Development 

The EDOP radar made it's first successful flight during the CAMEX experiments at the 
Wallops Flight Facility in September 1993. In this flight, the real-time reflectivity 
algorithm was tested, and the inital results are encouraging (Figure 5-1). Though the 
real-time Doppler algorithm was not tested in flight, it has been thoroughly simulated in 
the lab and appears to be properly implemented. Additional flights are planned at Wallops 
for the summer of 1994 in which the Doppler algorithms will be tested for the first time in 
the operational mode. 

Additional work must be completed before EDOP is fully meets the operational design 
goals specified in the original proposal. Perhaps the most important remaining task in 
completion of an automatic gain control (AGC) algorithm to increase the effective 
dynamic range of the system. This will be particularly important for Doppler processing 
because the linear coherent detectors are easily saturated. 

5.1 Automatic Gain Control 

One solution to this problem involves the implementation of hardware to perform 
sensitivity-time control and software to provide global gain control to compensate for 
gradual variations in the mean reflective echo. 

Sensitivity time control (STC) is a method used to avoid saturation of the receiver by 
strong returns at short ranges without compromising receiver sensitivity at longer ranges 
(Figure 5-2). After each pulse is transmitted, the receiver gain, which was initially greatly 
reduced, is increased with time to match the decrease in amplitude with range due to 
spreading. Over a range of 30 km, this amounts to an increase in receiver sensitivity by 
roughly 70 dB. Hence, maximum sensitivity is provided to enable detection of the 
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Convective clouds with a trailing stratiform region observed by the ER-2 
Doppler Radar (EDOP) on 5 October 1993 during the Convection and 
Atmospheric Moisture Experiment (CAMEX). The top panel shows 
reflectivity data (dBZ) from the forward looking antenna (30° off nadir) 
and the bottom panel shows the corresponding reflectivity data from the 
~ nadir viewing antenna. The forward image has been transformed to the 
coordinate system of the nadir looking antenna. 

The surface return is near 20km range and the ’’bright band” is near 15km. 

_ Each x-axis division is approximately 9km along the flight path. The thin 
line of enhanced reflectivity seen at 17km range on the forward image is 
an artifact that will be removed in future flights by modification of the 
radar hardware. 

Figure 5-1: Convective clouds observed during CAMEX (courtesy Jeff Caylor 



SENSITIVITY TIME CONTROL 


At low PRFs, saturation of the receiving system by strong 
return from short ranges is commonly avoided without loss of 
detection sensitivity at greater ranges through a feature called 
—sensitivity time control, STC. 



_ After each pulse has been transmitted, the system gain, which 
initially is greatly reduced, is increased with time to match the 
decrease in amplitude of the radar return with range. Maximum 
3am is usually reached well before the end of the interpulse 
""period. 

MAXIMUM RECEIVER GAIN 



Thereafter, the increase in sensitivity is continued by lowering 
the detection threshold until the noise limit is reached — i.e., to 
the point where the threshold is just far enough above the mean 
noise level to limit the false-alarm probability to an acceptable 
value. 



Thus, maximum sensitivity is provided at long ranges, where it 
is needed to detect the weak echoes of distant targets, while 
the strong return from short ranges is prevented from saturating 
the system. 

STC may be applied at various points in a system. From 
whatever point it is applied, it helps prevent saturation in all 
following stages. 




STC applied helps 

HERE HERE 


Figure 5-2: STC concept (Stimson, 1983) 






weakest echoes near the surface without saturating the receiver by strong echoes at higher 
altitudes. 

If necessary, the dynamic range of the system can be further increased by implementing 
digital automatic gain control (DAGC). This feature monitors the output of the A/D 
converters to build a continuously updated profile of the average receiver output. On the 
basis of the profile, embedded software produces a gain-control signal that is applied to a 
set of IF-attenuators ahead of the EF linear-detectors. A Kalman filter might be used to 
track the dc-level of the linear-detector output and center it within the dynamic range of 
the A/D converters. 

Additional resolution in range is possible by increasing the A/D sampling rate to 4 MHz 
from 2 MHz and narrowing the pulse width to 0.25 jis. These changes would allow an 
improved resolution of 37.5 meters or could be used to increase the number of effectively 
independent samples through range averaging. 

Additional test flights must be taken to thoroughly calibrate the radar system. These 
flights can be easily conducted up and down the coast of California over the ocean. 

5.2 Antenna Stabilization 

The Doppler velocity measured by an airborne radar is given by the dot product of the 
aircraft-velocity vector with the antenna-pointing vector. This product may then be 
related to the three-dimensional air-motion vector by transforming the aircraft body-fixed 
coordinate references to earth-fixed coordinates. Figure 5-3 shows the relation between 
the two coordinate systems. Heymsfield (1989) develops this transformation by first 
expressing the aircraft body-fixed coordinates as a linear transformation of the earth-fixed 
coordinates. It should be understood by the reader that such a transformation neglects the 
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2nd-order effects that result from the transformation between rotating and non-rotating 
reference frames. This section will instead assume that the aircraft undergoes rotation 
before the measurements are made. 


To describe the orientation of the aircraft relative to an earth-fixed coordinate system, it is 
sufficient to specify a three-dimensional transformation matrix that relates the two 
systems. The transformation is earned out by three consecutive rotations about the yaw-, 
pitch-, and roll-axis, respectively. Emphasis is placed on the term consecutive because 
the order in which they are carried must be consistent; these rotations do not obey a 
commutative property as demonstrated through Figure 5-4. 

Let the earth-fixed coordinate system be defined as the vector (x e , y e , z e ) and the aircraft 
body-fixed coordinate system be defined as the vector (x a , y a , z a ). 

Then it is possible to specify a rotation matrix A such that: 





y a 

= A 

Yc 

^ Z J 


v z J 


To transform the earth-fixed coordinates to the aircraft, the following rotations are applied 
in order: 

1) Rotate coordinate system (xl, yl, zl) about the z l axis over an angle vj/ as 
demonstrated in Figure 5-5. This angle is referred to as the heading or yaw angle. Note 
that earth-fixed system has been renamed from (x e , y e , z e ) to (x lf y lf Zj) for the purpose of 
this rotation. 
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X 2 = xi cos \|/ + y, sin vj/ 
y 2 = y, cos y - X i sin y 


X2 


’ cost)/ 

sin ij/~ 

Xi 

y 2 . 


—sin Vj/ 

COSVJ/ 

_y. 


2) Rotate coordinate system (x 2 , y 2 , z 2 ) about the y 2 axis over an angle 0 as demonstrated 
in Figure 5-5. This angle is referred to as the pitch angle. 

X 3 = X2cos0-Z2sin0 
Z3 = z2cos0 + X2sin0 


V 


COS0 

-sin0Tx 2 " 

- Z 3. 


sin0 

COS0 z 2 


3) Rotate coordinate system (x 3 , y 3 , z 3 ) about the x 3 axis over an angle <|> as demonstrated 
in Figure 5-5. This angle is referred to as the bank or roll angle. 

y = y 3 cos<t> + z3sin<l) 
z = -y 3 sint)) + Z3 cos<}> 


"y" 


" cost}) 

sin <}) If y 3 " 

z 


sin 4> 

cost)) z 3 


Therefore, the combined operation on the three-dimensional earth-fixed system appears as 


V 


"1 

0 

0 

COS0 

0 

-sin0 

COSVJ/ 

sin y 

0‘ 


y. 

= 

0 

COS()> 

sin <|> 

0 

1 

0 

-sinvy 

COS\|/ 

0 

y e 

.Z,_ 


0 

—sin <() 

COS<{) 

sin0 

0 

COS0 

0 

0 

1 

,Z e . 
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Figure 5 - 3 : Relationship between rotating and inertial 

coordinate systems (Roskam, 1979) 








For meteorological purposes, it is convenient to define the earth-fixed coordinate system 
such that x e points east, y e points north, and z e point up. Likewise it is convenient to 
define the aircraft body-fixed system such that x a points out the nose of the aircraft, y a 
points out the right wing, and z a points downward from the aircraft belly as depicted in 
Figure 5-3. Furthermore, the aircraft inertial navigation system reports the heading angle 
relative to North. The result of this is that, for zero angles of heading, pitch, and roll, the 
aircraft will be heading north. Therefore, the rotation matrices must be further modified to 
correctly transform the earth-fixed system to the body-fixed system. 

To do this, it is necessary to rotate the earth-fixed system by -90° in addition to the 
heading angle, y. This is accomplished by replacing \j / in the above matrices by -Vj/-90°. 
Next, one must pitch the earth-fixed system by 180° in addition to the pitch of the aircraft 
so that the z-axes of both systems will be correctly aligned. This is accomplished by 
replacing 9 in the above matrices by 0- 1 80°. These matrix modifications effectively align 
the two coordinate systems in the default state (no heading, pitch, or roll angles). 


X a 


"1 0 

0 

-COS0 

0 

sin0 

-siny 

-cosy 

0" 


y a 

= 

0 cos<}> 

sin<}> 

0 

1 

0 

COSVJ/ 

-siny 

0 

y e 

A- 


0 — sin (f) 

COS(J) 

-sin0 

0 

-COS0 

0 

0 

1 

_Z e _ 


The rotations about each axis are then combined to give the rotation matrix. A: 


A = 


cos 0 sin y 

sin 0 sin <J) sin y + cos <j> cos y 
cos 4> sin 0 sin y - sin <(> cos y 


COS 0 COS y 

sin 0 sin <)> cos y - sin y cos 0 
cos 4> sin 9 cos y + sin 4> sin y 


sin0 

-sin<|)cos0 

-cos<|>cos0 
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rotation SEQUENCE | 


ROTATION SEQUENCr ? 








The aircraft-motion vector can be expressed relative to the earth frame of reference by 


V ac — GSsin Tx + GScos Ty + w z 

c J e p e 


while the precipitation motion is described by 


V p = u, x e + v a y e + (w a + v t )z e 


where 

T 

aircraft track angle 

GS 

aircraft ground speed 

U a- V a 

horizontal winds at observation altitude 

W a 

air vertical velocity at observation altitude 

W P 

aircraft vertical velocity 

V « 

hydrometeor fallspeed at observation altitude 


Here, it is assumed that the hydrometeors follow the three-dimensional air motion and fall 

at their terminal velocities (v t <0). The radial velocity between the radar and precipitation 
is then given by 


v .=(V,-V„)e r 


For the case of a nadir-pointing antenna, e r = z a , reducing the Doppler velocity to 

v r = - w a cos P cos R - GS(sin P cos D + cos P sin R sin D) + 

+ w p cos P cos R - v, cos P cos R + ( u, sin H + v, cos H) sin P + (u 4 cos H - v, sin H) cos P sin R 
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Clearly, these terms result from the various vertical and horizontal air and aircraft motions. 
The effect of each term on the measured Doppler velocity is best illustrated by labeling 
and considering each term separately 


Vertical Air Motion: - w a cos P cos R 


Horizontal Aircraft Motion: 
Vertical Aircraft Motion: 
Fallspeed: 

Horizontal Air Motion: 


- GS(sin P cos D + cos P sin R sin D) 

+ w p cos P cos R 

- v, cos P cos R 

+ (u a sin H + v a cos H) sin P +(u a cos H - v a sin H) cos Psin R 


These variables can be determined both directly and indirectly from the information 
provided by the aircraft Inertial Navigation System. Consider first, a case whereby the 
antenna is fixed rigidly to the aircraft frame. It is desired to extract w a from the Doppler 
velocity. Though the Doppler velocity will be contaminated by terms related to the 
horizontal and vertical aircraft motion, these terms can be accounted for using information 
from the INS. The hydrometeor fallspeed can additionally be removed by estimating the 
fallspeed from an empirical relation of radar reflectivity to fallspeed. However, the motion 
of the horizontal winds is unknown for any observed altitude below the aircraft. 

Therefore, there is an error in the estimated w a which can be significant for strong 
horizontal winds. 


If the antenna is instead stabilized at nadir, these errors can be removed. In this case. 


P«0, 
R = 0 


and the measured Doppler velocity reduces to 
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V r =-( W a + V,) + W p 


Again, the hydrometeor fallspeed can be removed using an empirical relation. If the 
aircraft vertical speed, w p , can be obtained from the inertial navigation system, the mean 
precipitation vertical velocity can be extracted from the measured Doppler velocity. 


The inertial navigation system on the U-2 provides angle measurements that are accurate 
to within 15 seconds of arc or l A°. Neglecting the uncertainties in the radial velocity 
estimate as described in §4.2.3, the uncertain angle measurements affect the measured 
Doppler velocity in the following way: 


Table 5-1: Typical Flight Parameters 


Variable 

IEM5S2S9 

EEBBW 

P 

r 

0.06° 

R 

i° 

0.06° 

H 

5° 

0.1° 

T 

10° 

O 

o 

D 

5° 

0.3° 

GS 

200° 

2.0 ms' 1 

w„ 

0 ms' 1 

0.5 ms' 1 



0 ms -1 




20 ms 1 




0 ms -1 



Table 5-2: Motion Contributions to Doppler Velocity 


Doppler Velocity Term 

Calculated Contribution 

Calculated Uncertainty 

Vertical Air Motion 

5.0 ms -1 

Oms -1 

Horizontal Aircraft Motion 

2.28 ms* 1 

0.12 ms-‘ 

Vertical Aircraft Motion 

0 

0.5 ms-‘ 

Fallspeed 



Horizontal Air Motion 1 

0.03 ms-‘ 

0.0 ms -1 

Horizontal Air Motion 2 

0.35 ms -1 

0.0 ms- 1 

Total 

7.66 ms -1 

0.62 ms-' 
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Here, the fallspeed contribution to Doppler velocity and its uncertainty have been 
neglected. If the aircraft-motion contributions are subtracted from the measured Doppler 
velocity, the resulting vertical air-motion estimate is 5.38 ms 1 . It should be evident that 
the horizontal air-motion term is the significant contributor to this bias. As stated 
previously, because the horizontal winds at the observation altitude are almost never 

known, the only accurate way to remove this bias is through stabilization of the nadir 
antenna. 
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Appendix A - RIB Code 
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ribcwlox.s 


/* 

* File 

* 

* Purpose : This program acquires temperature information 

* each pulse and sends it to the VME queues. 

* After every 1000 pulses, the transmitter is 



* 

* 

switched to calibrate mode. 

— 

* Calls : 

* 




* Author 

: Shaun R. Nicholson 


— 

* 

Center For Research, Inc. 


* 

University of Kansas 



* 

* 

(913) 864-4835 



* Registers: rl => 0x610000 

; Board Status 


* 

r2 => 0x630000 

; Radar Status 


* 

r3 => counter 

; calibration counter 


* 

r4 => mi sc 

; Status register addr's 


* 

r5 => Attenuate 

; Calibration Attenuator tal 


* 

r6 => misc 



* 

rl => 



* 

r8 => 



* 

A 

II 



* 

rlO => BCR 

; Board Control Register 

— 

* 

rl 1 => 0x264000 

; VME A output queue 


* 

r 1 2 => Num_Pulses 

; # pulses/dwell - 2 


* 

rl3 => Dwell_cnt 

; Current Dwell # 

— 

* 

rl4 => One 

; 1.0 


* 

rl 5 => misc 

; Pulse # 


* 

rl6=> timer 

; for standby status 


* 

rl7 => timer 

; for standby status 


* 

rl8 => misc 



* 

r 19 => misc 

; Attenuator table offset 


* 

r20 => misc 

; misc status register 


* 

r21 => 



* 

r22 => Service 

; Interrupt Service Routine 


* 

* 

* 


* Revision 

* History : 

* 

* (0.x = Beta) 

* Revision: By: Date: Description: 


83 



S.R. Nicholson 07-28-93 Original. 


* 0.0 
* 

*/ 


#include "d:\shaun\dsp32c\bin\zone_addr.h" 

#include "d:\shaun\dsp32c\bin\rib_addr.h" 

.global main, Start, Serv_chk,Wait_rad,waitJoop, Wait Jnt,Service 
.global done,Dwell_cnt,One,temp,Num_Pulses,Cal_flag,Attenuate 

.rsect ".text" 
main: 

/* 

* (1) Disable Interrupt 1 (Pulse transmitted) 

* (2) Disable DMA transfers 

* (3) Set external memory partition A to 0 wait states 

* (4) Set external memory partition B to 2+ wait states 
*/ 

rl = OxOOOd /* No interrupts till trigger starts */ 
nop 

pew = r 1 /* interrupts and memory wait states */ 

ioc = 0x0 /* DMA */ 

r22e = Service /* Interrupt service routine */ 
goto Start 
nop 

.rsect ".R0" 

/* Initialize pointers */ 

Start: rle = BSR 

r2e = RSR 
r3e = 660001 
rile = 0x264000 
rl2e = Num_Pulses 
rl3e = Dwell_cnt 
rl4e = One 
rl 5 = *rl2 
rl6 = 32018 
rl7 = 937 


/* Board status register */ 

/* Radar status register */ 

/* #pulses + 1 BETWEEN calibrate cycles */ 
/* VME queue B */ 

/* Dwell # */ 

/* Start at Dwell #1 */ 

/* # pulses per dwell - 2 */ 

/* standby status timer */ 

/* timer */ 


I* 

* 

* Now, wait for the pilot to flip the operation switch from 

* STANDBY to RADIATE. While in STANDBY mode, status will be 

* recorded approximately once every 30 sec. 
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/* read RSR */ 


* 


*/ 

Wait_rad: rl8 = *r2 


nop 

r 1 8 = *r2 
nop 

r6 = rl8 


/* must read twice */ 


r6 = r6 & 0x0200 /* check STDBY_RAD bit */ 

if(ne) pcgoto Begin_trig /* bit lo - standby */ 

nop 

if(rl6— >=0) pcgoto Wait_rad /* millisecond timer */ 
nop 

if(r 1 7— >=0) pcgoto Wait_rad /* time for status? */ 

r 16 = 32018 /* regardless, reload inner timer */ 

call Status (r8) /* collect status */ 

nop 

r 1 7 = 937 /* reload outer timer */ 

goto Wait_rad 

nop 


/* 

* The pilot has flipped the switch to the RADIATE position. 

* The trigger will now be turned on... 

* 

*/ 


Begin_trig: rlOe = BCR 

r6 = 0x00a3 
*rl0 = r6 
r6 = 0x80a3 
*rl0 = r6 
r6 = 0x00a3 
*rl0 = r6 
rl8 = 0x800d 
pew = r 1 8 

/* 


/* Board Control register */ 

/* the radiate signal */ 

/* has been captured - */ 

/* start trigger now */ 

/* trigger started */ 

/* Now start interrupts! */ 

/* interrupts and memory wait states */ 


Once the trigger is provided, the A/D's will begin dumping samples 

* into their queues. The E-6 will in turn start dumping results 

* to the hard drive for storage. At this point, the program simply 

* waits for interrupts and checks the fault registers for an 
indication of trouble or transmitter turnoff. An interrupt will 

* occur with the transmission of each pulse - at the end of each 
dwell, sensor and status information are sent to the VME queue B 

* for storage on physical media. 

*/ 
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/* read RSR */- 


Wait_int: rl8 = *r2 
nop 

rl8 = *r2 
nop 

r6 = r 1 8 

r6 = r6 & 0x0200 
if(ne) pcgoto Wait Jnt 
nop 

call Timer (r8) 
nop 

r 1 8 = OxOOOd 
pew = r 1 8 

rl8 = *r2 
nop 

rl8 = *r2 
nop 

r6 = rl8 

r6 = r6 & 0x0200 

if(eq) pcgoto terminate 
nop 

rl8 = 0x800d 
pew = r!8 

goto Waitjnt 
nop 


/* must read twice */ 

/* check STDBY_RAD bit */ 
/* still transmitting... */ 

/* wait a bit */ 

/* stop interrupts */ 

/* read RSR */ 

/* must read twice */ 

/* check STDBY_RAD bit */ 
/* transmitting */ 

/* restart interrupts */ 


/* If FAULT2 (hard fault) was asserted, the transmitter has failed during 

* both the main start and attempted restart. In this case, (or if 

* the transmitter was turned off by the pilot), the GO/STOP bit should 

* be set to one thereby indicating that the E-6 should close its 

* data files and terminate. 

*/ 

terminate: rl8 = 0x4003 /* Signal E-6 to terminate *! 

*rl0 = rl8 

NextLife: pcgoto NextLife 
nop 

/* SUBROUTINE: STATUS 
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* Description: This subroutine is used to acquire status information 
(temperature and pressure) about the radar transmitter and receiver. 
These data are collected periodicly and sent to the hard drive via 

* VME bus B. 

*/ 

Status: rl4e = temp 

r20 = *r 1 /* Board Status register */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 

r20 = *r2 /* Radar Status register */ 

nop 

r20 = *r2 /* must read twice */ 

nop 

*rl4 = r20 


aO = (*rl l=*rl4) + aO 


r4e = ADCDACO 
r20 = *r4 

/* temp 1 */ 

nop 

r20 = *r4 

/* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 


r4e = ADCDAC1 
r20 = *r4 

/* temp 2 */ 

nop 

r20 = *r4 

/* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 


r4e = ADCDAC2 
r20 = *r4 

/* temp 3 */ 

nop 

r20 = *r4 

/* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 


r4e = ADCDAC3 
r20 = *r4 

/* temp 4 */ 
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/* must read twice *! 


nop 

r20 = *r4 
nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 

r4e = ADCDAC4 
r20 = *r4 /* temp 5 */ 

nop 

r20 = *r4 /* mus t re ad twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
r4e = ADCDAC5 

r20 = *r4 /* temp 6 */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 

r4e = ADCDAC6 
r20 = *r4 
nop 

r20 = *r4 
nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 

r4e = ADCDAC7 
r20 = *r4 
nop 

r20 = *r4 
nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
r4e = ADCDAC8 

r20 = *r4 /* rib temp */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*r!4) + aO 


/* pressure 2 */ 

/* must read twice */ 


/* pressure 1 */ 

/* must read twice */ 
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r4e = One 

aO = (*rl l=*rl3) + aO /* Dwell # */ 

return (r8) 
nop 


/* SUBROUTINE: TIMER 
* 

* Description: This subroutine provides a 10 second time delay. 
*/ 


Timer: r6e = 10000000 /* about 10 seconds */ 

/* ( 1 microsecond/loop) */ 
waitjoop: r6e = r6 - 1 /* tick lock, tick tock... */ 

nop 
nop 
nop 
nop 
nop 
nop 
nop 

if(ne) pcgoto waitjoop /* done? */ 
nop 

return (r8) 
nop 


/* 

* This is the interrupt service routine. When an interrupt occurs, 

* the pulse # counter is decremented. Once it goes to zero (ie. 

* at the end of each dwell), status information is transferred to 
the VME queue for storage and the counter reset 

*/ 

Service: r20 = OxOOOd /* Mask interrupts */ 

nop 

r4e = Cal_chk 
pew = r20 

if(rl5~>=0) goto r4 /* no status info this time */ 
nop 
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/* Found last pulse this dwell - record status information */ 
Serv_chk: rl4e = temp 


— . 

r20 = *rl 
nop 

*rl4 = r20 

/* Board Status register */ 

— 

aO = (*rl l=*rl4) + aO 



r20 = *r2 

/* Radar Status register */ 

— 

nop 



r20 = *r2 
nop 

/* must read twice */ 


*rl4 = r20 

aO = (*rl l=*rl4) + aO 


— 

r4e = ADCDACO 



r20 = *r4 
nop 

/* temp 1 */ 


r20 = *r4 
nop 

*rl4 = r20 

/* must read twice */ 


aO = (*rl l=*rl4) + aO 



r4e = ADCDAC1 



r20 = *r4 
nop 

/* temp 2 */ 

— ■ 

r20 = *r4 
nop 

*rl4 = r20 

/* must read twice */ 


— 

aO = (*rl l=*rl4) + aO 



r4e = ADCDAC2 



r20 = *r4 

/* temp 3 */ 


nop 

r20 = *r4 

/* must read twice */ 


nop 

*rl4 = r20 




a0 = (*rll=*rl4) + a0 



r4e = ADCDAC3 



r20 = *r4 

/* temp 4 */ 


nop 

r20 = *r4 

/* must read twice */ 


nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
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r4e = ADCDAC4 

r20 = *r4 /* temp 5 */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
r4e = ADCDAC5 

r20 = *r4 /* temp 6 */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
r4e = ADCDAC6 

r20 = *r4 /* pressure 1 */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
r4e = ADCDAC7 

r20 = *r4 /* pressure 2 */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 
r4e = ADCDAC8 

r20 = *r4 /* rib temp */ 

nop 

r20 = *r4 /* must read twice */ 

nop 

*rl4 = r20 

aO = (*rl l=*rl4) + aO 

aO = (*rl l=*r!3) + aO /* Dwell # */ 

rl4e = *rl3 

nop 

r!4e = r 14 + 1 /* Increment Dwell # */ 
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*rl3 = rl4e 
rl5 = *rl2 

Cal_chk: r4e = Cal_flag 

r20 = *r4 
nop 

if(ne) pcgoto Calibrate 
nop 


/* reset pulse counter */ 

/* are we now calibrating? */ 

/* yup! */ 


r3e = r3 - 1 
if(ne) pcgoto done 
nop 

Cal_Start: r4e = RCRO 

r20 = 0x16 
*r4 = r20 
rl9 = 0 
r4e = Cal_flag 
*r4 = r20 

Calibrate: r5e = Attenuate 

r4e = RCR 1 
r5e = r5 + rl9 
r20 = *r5 
nop 

if(eq) pcgoto nomore 
nop 

*r4 = r20 
rl9 = r 19 + 2 


/* time to begin calibration? */ 
/* not yet! */ 


/* Radar Control reg. */ 

/* 16=2200 Hz, 12=4400 Hz */ 

/* turn on calibrate mode */ 

/* offset into table */ 

/* set flag */ 

/* Table base address */ 

/* Radar Control reg. */ 

/* add offset */ 

/* get desired setting */ 

/* end of table? */ 

/* nope, set attenuators! */ 

/* point to next entry */ 


/* finished with current calibration dwell */ 
cal_done: goto done 
nop 


/* finished with all calibration dwells, return to normal processing */ 
nomore: r4e = Cal_flag 


*r4 = r20 
r4e = RCRO 
r20 = Oxle 
*r4 = r20 
r3e = 660001 


/* clear calibrate flag */ 

/* Radar Control reg. */ 

/* end calibrate mode */ 

/* le=2200 Hz, la=4400 Hz +/ 

/* #pulses + 1 BETWEEN cal cycles */ 


done: 


if(ireq 1 Jo) pcgoto done /* wait till irq returns */ 
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nop 

r20 = 0x800d 

pew = r20 

rl3e = Dwell_cnt 

rl4e = One 

ireturn 

nop 


/* high before exiting */ 
/* restore interrupts */ 

/* restore pointer */ 

/* restore pointer */ 


.rsect ".text" 
.align 4 

Dwell_cnt: 

fltbits 0x0 

/* current Dwell # */ 

One: 

float 1.0 


temp: 

fltbits 0x0 


Num_Pulses: 

int 6598 

/* # xmit pulses/status - 2*1 

Cal_flag: 

int 0x0 

/* Calibration dwell flag */ 

Junk_variable: int Oxface 

.align 4 

Attenuate: 

220*int OxOfcO 



220* int OxOfcl 
220*int 0x0fc2 
220*int 0x0fc4 
220*int 0x0fc8 
220*int OxOfdO 
220* int OxOfeO 
220*int OxOfcO 
int 0x0000 

/* end of table */ 
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Appendix B - ACQ Code 
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/* 

* File acqreal.s 

* Purpose : This is the data scattering program to be flown 

during the Wallops Island experiments during 
September 1993. For each transmitted pulse, 

436 tags are supplied to route the raw data to 
LUA Board 1 for reflectivity processing. It 
is intended that the system operate at a PRF of 
2200 Hz, 0.25 usee pulse width, and 1 MHz sampling. 


* Author : S.R. Nicholson 

* Date : 8/10/93 

*/ 


#include "d:\shaun\dsp32c\bin\acq_addr.h" 

.rsect ".text" 

.align 4 
main: 

/* 

* (1) Enable interrupt one from dataqaffn 

* (2) Disable DMA transfers 

* (2) Set external memory partition A to 0 wait states 

* (3) Set external memory partition B to 2+ wait states 
*/ 


ioc = 0x0 
goto Start 
nop 

.rsect ”.R0" 

.align 4 

Start: rle = ADDQAB 

r3e = Range_Gates 
r5e = ADCTL0 
r6 = OxOOOd 
r7 = 0x800d 
rl le = Int_flag 
rl3 = 0 

rl4 = 0x0 
rl6 = 998 
r21 = 1 
r22e = Intsvc 
pew = r7 


/* DMA */ 


/* VME output queue's */ 

/* # range_gates */ 

/* to control A/D's */ 

/* mask Int_l */ 

/* enable Int_l */ 

/* A/D control bit */ 

/* 0 = ARM, 1 = Stop */ 

/* Always! */ 

/* #good - 2 before token pass */ 

/* Always! */ 

/* load address of qempty sve routine */ 


C.-1, 


95 



Send_Good: call Write_Tags (r8) 
nop 


Monitor: r20 - *rl 1 /* W ait for interrupt *1 

nop 

if(eq) pcgoto Monitor 
nop 

*rl 1 = r 14 /* clear interrupt flag */ 

if(r 16 — >=0) pcgoto Send_Good /* time to pass token? */ 
nop 


r 1 6 = 998 

/* reload pulse counter */ 

r8e = First_Time 

/* don't pass the first time 

r20 = *r8 


nop 


if(eq) pcgoto Send Good 


*r8 = r21 

/* set flag regardless */ 

Pass_Token: r2e = Pass_Cmd 


aO = (*rl=*r2) + aO 

/* pass token to LUA */ 


goto Send_Good 
nop 


/* subroutine 

- write valid tags for 1 pulse */ 


Write_Tags: 

r 1 8 = *r3 

/* #ranges */ 

Writejoop: 

r2e = Keep_Tag 
do 0,rl8 



aO = (*rl=*r2++) + aO 

/* Tag 1 */ 


nop 

*rl 1 = rl4 

/* clear Int_flag */ 


return (r8) 
nop 


/* Service routine for external interrupt one: this interrupt is 

* generated each time the pre-trig goes low. During the interrupt, 

* the Trigger Source bit in ADCTLO is pulled low, thereby arming 

* the A/D's to begin conversions once the pre-trig returns high. 

* The DSP then burst writes destination tags for each range gate 

* to be acquired, increments a dwell counter, and reloads a 

* range-gate count down counter that is used in interrupt 6. 
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*/ 

Intsvc: pew = r6 /* mask Interrupt 1 */ 

*r5 = r 1 3 /* reload A/D control bit */ 

*rl 1 = r21 /* set Int_flag */ 

Exit_l : if(ireq ljo) pegoto Exit_l 

nop 

pew = r7 /* re-enable interrupt 1 */ 

iretum 


.rsect ".text" 

.align 4 

Range_Gates: int 217 
Int_flag: int 0x0 

First_Time: int 0x0 


/* # ranges - 1 */ 

/* set high in Int. 1 */ 
/* first dwell flag */ 


.rsect ".Rl" 

/* Format : Data_B : Data_A : Dest_B : Dest_A: */ 
.align 4 

Pass_Cmd: fltbits 0x0000 1 f 1 f 

Keep_Tag: 10*fltbits 0x00000010 

10*fltbits 0x00000011 
10*fltbits 0x00000012 
10*fltbits 0x00000013 
200* fltbits 0x00000000 


97 



Appendix C - Reflectivity Code 
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/* 


* File 

* 

refzoneO.s 

* Purpose : 

This program is used for speed testing of the 

* 

reflectivity algorithm. It is nearly identical 

* 

to reflect, s. 

* Calls : 

* 


* Author : 

Shaun R. Nicholson 

* 

Center For Research, Inc. 

* 

University of Kansas 

* 

(913) 864-4835 


* 

* Registers: rl 

* r2 

* r3 

* r4 

* r5 

* r6 

* r7 

* r8 

* r9 

* rlO 

* rl 1 

* rl2 

* rl3 

* rl4 

* rl 5 

* rl 6 

* rl7 

* r 1 8 

* rl9 

* r20 

* r21 

* r22 

* 

* 

* 

* Revision 

* History : 

* 

* 

* 

* 

* 


=> QDATA 
=> QDATB 
=> DESTA 
=> DESTB 
=> NadirXMean 
=> NadirCoMean 
=> FwdCoMean 
=> Range_Gates 
=> 

=> Nadir_Hold 
=> Fwd_Hold 
=> misc 
=> misc 
=> misc 
=> misc 
=> misc 
=> misc 

=> Range count 
=> 

=> Pulse count 
=> 

=> 


(0.x = Beta) 

Revision: 

0.0 


; DSP Queue A 
; DSP Queue B 
; Bus A Destination Register 
; Bus B Destination Register 
; Nadir X-Pol Echo Power Means 
; Nadir Co-Pol Echo Power Means 
; Forward Co-Pol Echo Power Means 
; Addr containing # range gates 

; Input Buffer for Nadir Channel 
; Input Buffer for Forward Channel 


; Used to output results 
; Used for results dest. tags 


; # range gates 
; # pulses/dwell 


By: 

S.R. Nicholson 


Date: 

07-24-92 


Description: 

Original. 
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*/ 

/* #include "d:\shaun\dsp32c\bin\zone_addr.h" */ 


.global main, Start, Next_Dwell,Next_Block,Next_Pulse,Compress,Output 
•global NadirXMean,NadirCoMean,FwdCoMean 
•global Count2,Num_Block,Num_Pulses,Range_Gates,Comp_Factor 
•global Nadir_Hold,Fwd_Hold,PassToken 
•global QDAT A,QDATB ,DESTA,DESTB 

.rsect ".text" 
main: 

/* 

* ( 1 ) Disable interrupts 

* (2) Disable DMA transfers 

* (2) Set external memory partition A to 0 wait states 

* (3) Set external memory partition B to 2+ wait states 
*/ 

rl = OxOOOd 

pew = rl /* interrupts and memory wait states */ 

ioc = 0x0 /* DMA */ 

goto Start 
nop 


.rsect ”.R0" 

/* Initialize pointers */ 

Start: rle = QDATA 

r2e = QDATB 
r3e = DESTA 
r4e = DESTB 
r8e = Range_Gates 
Next_Dwell: rl4e = Dwell_cnt 

rl 5 = *rl4 


/* Bus A DSP Queue */ 

/* Bus B DSP Queue */ 

/* Bus A Destination Register */ 

/* Bus B Destination Register */ 

/* Addr containing # range gates */ 
/* Update Dwell count */ 


nop 

r 1 5 = rl5 + 1 
*rl4 = rl5 
rl4e = Range_cnt 
rl5 = 0 
*rl4 = rl5 
r5e = NadirXMean 
r6e = NadirCoMean 
r7e = FwdCoMean 


/* increment count */ 
/* and save... */ 

/* next, zero out */ 

/* the range count */ 


* Clear arrays used for holding partial sums and resultant means 
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* for each range gate. 

*/ 

rl3e = Zero 
rl8 = *r8 
nop 

do 2,rl8 

aO = (*r5++=*rl3) + aO 
aO = (*r6++=*rl3) + aO 
a0 = (*r7++=*rl3) + a0 


/* N range gates */ 

/* Clear NadirXMean array */ 
/* Clear NadirCoMean array */ 
/* Clear FwdCoMean array */ 


/* 

* Now, read a pulse of reflectivity data from the input FIFOs and 

* store in separate buffers. The data for an entire dwell will be 

* read pulse by pulse. 

* 

* This program expects that each word read from FIFO A contains 

* a single range gate from the Nadir Cross-polarization receiver 

* in the upper 16 bits, and a single gate from the Nadir Co-polarization 

* receiver in the lower 16 bits. Each word read from FIFO B is 

* expected to contain a gate from the Forward Co-polarization 

* receiver in the lower 16 bits. 

* 

* For each range gate, the pulses in a given dwell are summed 

* and multiplied by a scaling factor (Comp_Factor). They are then 

* stored in the mean-reflectivity table. Data for 

* a dwell are simultaneously summed and scaled, taking full 

* advantage of the DSP32Cs multiply-accumulate architecture. 

*/ 


rl2e = Num_Pulses 
r20 = *rl2 
r5e = NadirXMean 
r6e = NadirCoMean 
r7e = FwdCoMean 
rl5 = 4 


/* # pulses/dwell */ 


/* offset */ 


Next_Pulse: r 1 8 = *r8 

rlOe = Nadir_Hold 
rl le = Fwd_Hold 
do l,rl8 

aO = (*rlO++=*rl) + aO 
aO = (*rl l++=*r2) + aO 


/* get # range gates */ 


/* read pulse from input FIFOs */ 
/* Read one range from Nadir */ 
/* and one from Forward ant. */ 


r 18 = *r8 /* get # range gates */ 

rlOe = Nadir_Hold 
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rile = Fwd_Hold 
do 5,rl8 

al = float(*rlO++) 
aO = float(*rlO++) 
a2 = float(*rl l++rl5) 
*r6++ = a3 = *r6 + al 
*r5++ = a3 = *r5 + aO 
*r7++ = a3 = *r7 + a2 
r5e = NadirXMean 
r6e = NadirCoMean 


/* Nadir Co */ 

/* Nadir X */ 

/* Forward Co */ 

/* measurements for all */ 

/* accumulate echo power */ 
/* three channels */ 

/* reset pointers */ 

/* for next pulse */ 


if(r20— >=0) pcgoto Next_Pulse 
r7e = FwdCoMean 

/* more pulses this dwell? */ 

/* (this one is executed too!) */ 

Compress: rl8 = *r8 

/* get # range gates */ 
/* Scaling factor */ 

r!2e = Comp_Factor 
do 2,rl8 

*r5++ = aO = *r5 * *rl2 
*r6++ = aO = *r6 * *rl2 
*r7++ = a0= *r7 * *rl2 



/* 

* Output results: 

* 


Results are output as follows: Nadir polarisations arc output 
via Bus A with the cross-polarisation in the upper 16 bits and 

* * he co ‘P olansatl °n in the lower 16 bits. Forward Co-polarisation 
is sent in the lower 16 bits of Bus B. Bits 16-25 (Bus B) contain 
^range gate # while bits 26-3 1 contain the Dwell # (modulo 64) 


Output: r5e = NadirXMean 

r6e = NadirCoMean 
r7e = FwdCoMean 
rl8 = *r8 
rl2e = pack 
rl3e = Range_cnt 
rl6 = 0x44 
*r3 = rl6 
*r4 = r 1 6 
do 26,rl8 

♦rl2++ = aO = int(*r6++) 
*rl2- = a0 = int(*r5++) 
r 14e = Dwell_cnt 


/* get # range gates */ 

/* pack X & Co polarisations */ 

/* tag current range */ 

/* send results to Board 4, zone 4 */ 
/* Bus A destination */ 

/* Bus B destination */ 

/* get Nadir Co */ 

/* get Nadir X */ 
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/* Get current dwell # */ 


rl5 = *rl4 
nop 

aO = (*rl=*rl2) + aO 
r 15 = rl5 & 0x3f 
rl5 = rl5»>l 
r 15 = rl5>» 1 
rl5 = rl5»>l 
r 15 = rl 5»> 1 
rl5 = rl5»>l 
rl5 = r 15»> 1 
rl5 = rl5>»l 
rl4 = *rl3 
nop 

rl4 = r 14 + 1 
*rl3 = rl4 
rl4 = rl4 & 0x3ff 
rl5 = r 1 5 I rl4 
*rl2++ = aO = int(*r7++) 
nop 

*rl2 — = r 15 

nop 

nop 

nop 

aO = (*r2=*rl2) + aO 
nop 

/* 

* Pass the token 
*/ 

PassToken: r 1 3 = Ox 1 f 

*r3 = r 1 3 
*r4 = r 1 3 
r 13 = 0x11 
*rl = rl 3 
*r2 = rl 3 
goto Next_Dwell 
nop 

.rsect ".text" 

.align 4 

NadirXMean: 436*float 0.0 

NadirCoMean: 436*float 0.0 

FwdCoMean: 436*float 0.0 

Nadir_Hold: 436*fltbits 0x0 


/* send Nadir polarisations */ 
/* truncate to 6 bits */ 

/* shift right 1 bit */ 

/* shift right 2 bits */ 

/* shift right 3 bits */ 

/* shift right 4 bits */ 

/* shift right 5 bits */ 

/* shift right 6 bits */ 

/* shift right 7 bits */ 

/* get previous range # */ 

/* update it */ 

/* and save it */ 

/* truncate to 10 bits */ 

/* and pack it in place */ 

/* Forward Co in lower 16 */ 

/* store in upper 16 bits */ 

/* latency... */ 

/* latency... */ 

/* latency... */ 

/* ship it! */ 

/* latency */ 


/* Pass token cmd */ 

/* Bus A */ 

/* Bus B */ 

/* pass to node 1 , zone 1 */ 


/* Mean Reflectivity table */ 
/* # range gates */ 


/* Holding area for packed reflectivity */ 
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Fwd_Hold: 


/* measurements from nadir antenna */ 

/* # range gates */ 

436*fltbits 0x0 /* Holding area for packed reflectivity */ 

/* measurements from forward antenna */ 
/* # range gates */ 


.rsect ".Rl" 


.align 4 


Zero: 

float 0.0 

Num_Pulses: 

int 1998 

Range_Gates: 

int 108 

.align 4 


pack: 

float 0.0 

holdpack: 

int 0 

PackFlag: 

int 0 


Fubar: 

int Oxface 

stopflag: 

int Oxface 

.align 4 

addrflag: 

float 0.0 

Comp_Factor: 

float 0.0005 

.align 4 

Range_cnt: 

int 0 

Dwell_cnt: 

int 0 

QDATA: 

float 0.0 

QDATB: 

float 0.0 

DESTA: 

float 0.0 

DESTB: 

float 0.0 


/* # pulses/dwell - 2 */ 

/* Number of range gates - 1 */ 


/* holding spot for result packing */ 

/* holds extra result between range packing */ 
/* >0 => SNR from last range waiting to be 

* packed. It is sitting in holdpack. 

* =0 => previous range completely packed */ 
/* For test purposes only */ 

/* For test purposes only */ 

/* for test purposes only */ 

/* l/(# pulses/dwell) */ 

/* for test purposes */ 

/* for test purposes */ 
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Appendix D - Autocovariance Code 
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File : ppzoneO.s 

Purpose : This program implements the pulse-pair algorithm. 

It is intended to be run as zonel on an LUA200 
designated as Board ID = 1. Results are packed. 

Calls : 

Author : Shaun R. Nicholson 

Center For Research, Inc. 

University of Kansas 
(913)864-4835 

Revision 
History : 


(0.x = Beta) 
Revision: By: 

Date: 

Description: 

0.0 

S.R. Nicholson 

11-15-91 

Original. 

0.1 

S.R. Nicholson 

01-09-92 

Vp calculations were 

0.2 

S.R. Nicholson 

03-10-92 

not done correctly. 
Token passing was not 

0.3 

S.R. Nicholson 

03-31-92 

done correctly. 

Removed FIFO allocation. 
Replace if() goto <addr> 

0.4 

S.R. Nicholson 

04-07-92 

with if() goto rxx for 
24-bit addressing, 
ver. 0.3 was corrupt 

0.5 

S.R. Nicholson 

06-22-92 

Replace if() goto rxx 
with if() pcgoto <addr> 
to use offset addressing. 
Reduce #pulses/range gates 

0.6 

S.R. Nicholson 

07-01-92 

for token pass debugging. 
Stop after one dwell 

0.7 

S.R. Nicholson 

07-02-92 

Count 1 defined wrong. 

0.8 

S.R. Nicholson 

07-12-92 

Remove code to save 


current DSP output 
queue. This was used 
only for simulation. 



.global main,Dwell,Block,Range,IQstart,RStart,Vpstart,Out 
.global _atan,_div,_loge,_xtoy 

.global QDATA,Rout,DESTA,IpPrev,Ip,QpPrev,Qp,NumPulse,RangeGates,Vpcoeff 
.global V arCoeff,OneThird,FourThirds,Zero,temp 1 ,temp2,temp3,temp4,Count 1 
.global Count2,Count3,Count4,Block,IQPack,I0,Il,I2,Ql,Q2,R0,Rl,R2,Vp,Var,SNR 
.global Routaddr,PartSum,theta,pack,holdpack,PackFlag,xtra,noxtra 
.global Fubar 

#include "c :\dosapps\dsp3 2c\bin\zone_addr.h " 

.rsect ".text" 
main: 

/* 

* (1) Disable interrupts 

* (2) Disable DMA transfers 

* (2) Set external memory partition A to 0 wait states 

* (3) Set external memory partition B to 2+ wait states 
*/ 

rl = OxOOOd 

pew = rl /* interrupts and memory wait states */ 

ioc = 0x0 /* DMA */ 

goto Dwell 
nop 

/* 

* At the beginning of each dwell, memory registers used to accumulate 


* results need to be cleared.. 
*/ 

.rsect ”.R2" 


Dwell: r2e = Dwell_cnt 

r 14 = *r2 
nop 

/* Update dwell count */ 

rl4 = rl4 + 1 

/* increment count */ 

*r2 = rl4 

/* and save... */ 

r2e = Range_cnt 

/* next, zero out */ 

r 14 = 0 
*r2 = rl4 

r2e = PartSum 
rl4e = Count4 
rl = *rl4 
r3e = Zero 
do 0,rl 

/* the range count */ 
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aO = (*r2++=*r3) + aO 

rl4e = NumBlock /*# blocks per dwell - 2 */ 
rl8 = *rl4 


/* 


* Read 32 bit data from the input FIFO and store in external memory. The 

* data are packed I and Q (12 bits each sign extended to 16 bits) with I 

* in the most significant 16 bits. These data will be transferred directly 

* to memory and unpacked later. This algorithm will process a variable 

* number of range gates. To allow this, the SBC must load the required number 

* of iterations to be done into the 16 bit Zone memory variable Count 1. 

* Count 1 = #range gates * #pulses each range - 2 
*/ 


Block: 


loop: 


/* 


rle = QDATA 
r2e = IQPack 
r!4e = Count 1 


r3 - *rl4 /* read data from input FIFO */ 

aO = (*r2++=*rl) + aO 
if(r3— >=0) pcgoto loop 
nop 

rle = RangeGates /* #range gates to process */ 
rl7 = *rl 


rle = Routaddr 
r2e = QDATA 

*r 1 = r2e /* Save address of Output FIFO */ 

r 1 3e = IQPack /* Base address for packed I & Q */ 

rl2e = PartSum /* Partial results of complex math */ 


* At this point, a block of packed I & Q data has been read from the input 

* FIFO and stored in external RAM. Next, they will be separated into 

* I and Q components, converted to the special DSP32C floating point format, 

* and stored in on-chip RAM. In-phase data will be stored in internal RAM 

* bank 0 while Quadrature data will be stored in internal RAM bank 1. Note 

* that only sufficient data for one range gate is expanded at a time. These 

* are then processed and intermediate results stored. They cycle repeats for 

* the next range gate. This is due to memory limitations. 

*/ 


/* 

* Separate I and Q data, and convert to float... 

*/ 

Range: rle = Ip 

r2e = Qp 

r3e = r 1 3 /* Get base addr of data for range X */ 

r!3e = rl3+4 /* Point to base for next range */ 
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rl4e = Count3 /* get pulse-to-pulse distance */ 
rl5 = *rl4 

rl4e = Count2 /* # pulses each (partial) range */ 

r5 = *rl4 

nop 

do 1 ,r5 /* enough data for 1 partial range */ 

*rl++ = aO = float(*r3++) 

*r2++ = aO = float(*r3++rl5) 

/* 

* Define Pointers... 

*/ 

rle = rl2 /* Points to IpPrev stored in PartSum 

* for the current range */ 

r6e = IpPrev 
r7e = IpPrev+4 
r8e = IpPrev+8 
r9e = QpPrev 
rlOe = QpPrev+4 
rl le = QpPrev+8 

aO = (*r6=*rl++) + aO /* Restore IpPrev */ 

aO = (*r7=*rl++) + aO /* Restore IpPrev+1 */ 

aO = (*r9=*rl++) + aO /* Restore QpPrev */ 

aO = (*rlO=*rl++) + aO /* Restore QpPrev+1 */ 

/* Now, rl points to Previous 10 */ 
r2e = r 1 2+20 /* Previous II*/ 

r3e = r 12+24 /* Previous 12 */ 

r4e = r 12+28 /* Previous Q1 */ 

r5e = r 12+32 /* Previous Q2 */ 

/***+**++*+++++++*+ ************************************************** 

* 

* Accumulate pulse pairs... 

* 

* 10 += Ip A 2 + Qp A 2 

* II += Iplp+l + QpQp+1 

* Q1 += -IpQp+1 + Qplp+l 

* 12 += IpIp+2 + QpQp+2 

* Q2 += -IpQp+2 + QpIp+2 

* 

IQstart: rl4e = Count2 /* How many pulses? */ 

r 16 = *rl4 
nop 

do 12,rl6 

aO = *r6++ /* IpPrev */ 

al = *r9++ /* QpPrev */ 
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no P /* latency... */ 

a2 = *rl + aO*aO /* a2 = 10 + Ip A 2 */ 

*rl = a2 = a2 + al * al /* 10 = Ip A 2 + Qp A 2 *1 

a2 = *r2 + aO * *r7 /* a2 = II + Iplp+l */ 

*r2 = a2 = a2 + al * *rlO /* II = a2 + QpQp+1 */ 

a2 = *r4 - aO * *rlO++ /* a2 = Q1 - IpQp+1 */ 

*r4 = a2 = a2 + al * *r7++ /* Q1 = a2 + Qplp+l */ 

a2 = *r3 + aO * *r8 /* a2 = 12 + IpIp+2 */ 

*r3 = a2 = a2 + al * *rl 1 /* 12 = a2 + QpQp+2 */ 

a2 = *r5 - aO * *rl 1++ /* a2 = Q2 - IpQp+2 */ 

*r5 = a2 = a2 + al * *r8++ /* Q2 = a2 + QpIp+2 */ 

/* 

* Save last two I & Q pulses - they will be needed for 

* the next data block! 

*/ 

aO = (*rl2++=*r6) + aO /* First, Ip */ 

aO = (*rl2++=*r7) + aO /* then Ip+1 */ 

aO = (*rl2++=*r9) + aO /* Qp */ 

aO = (*rl2++=*rl0) + aO /* and Qp+1 */ 

nop 

r 1 2e = r 1 2 + 20 /* skip past copies of I0-Q2 

* stored in Parts um */ 

I* 

* More range gates to process this Block ? 

*/ 


if(rl7 — >=0) pcgoto Range /* Yup! */ 

nop 

r!2e = PartSum /* Nope, reset partial sum */ 

/* pointer for the next block */ 


/* 

* Is there another block of data to process this Dwell 
*/ 

if(r 1 8 — >=0) pcgoto Block /* Yup! */ 

nop 

/* 

* If not, then the complex summation is finished and it is 

* time to compute the various parameters of interest. 

*1 

rl4e = RangeGates 
rl7 = *rl4 
RStart: rle = 10 

rl2e = rl2 + 16 /* point to 10 for current range */ 

d° 0.4 /* copy I & Q sums to named memory 

* variables for further 
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* processing */ 


aO = (*rl++=*rl2++) + aO 

/****************************************** j^***********************^ 

/* R(O) = IO/N */ 

y******* ********** d********* ********************************* ********/ 

chk2: nop 

nop 
nop 
rle=I0 
r2e=R0 

a0=(*r2=*rl) + aO 

/*******************************> ************************************y 

/* IR(l)l = sqrt(Il A 2 + Ql A 2)/N */ 

y********************************************************##,,,,,,,,,,,,,,,,,,,,,,,,,,^ 

rle = II 

r2e = Q 1 

r3e = tempi 

aO = *rl * *rl 

*r3 = aO = aO + *r2 * *r2 

call _sqr (rl4) 

nop 

.align 4 

int24 temp3 
int24 templ.Rl 

.align 4 

/********* ******* ^^^^^^^i******************#******#***#*******^ 

/* IR(2)I = sqrt( I2 A 2 + Q2 A 2 )/N */ 

y***** ************************************** *** 1(c ,| t , tt , (c ***** J|tl) , 1(cl) ,*, ) , Jk , ( ,, |<!4l , |1 , 4l y 

rle = 12 

r2e = Q2 

r3e = temp 1 

aO = *rl * *rl 

*r3 = aO = aO + *r2 * *r2 

call _sqr (rl4) 

nop 

.align 4 

int24 temp3 
int24 templ,R2 

.align 4 

I**************************************************#**#*************/ 

/* Vp = VpCoeff * arg[R(l)] */ 

/ * * ** * * ** * * * * * * * Dc * * * * * ** * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * i«e *y 
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Vpstart: 

rle=Il 



r2e=temp 1 
aO=-*rl 

- 


*r2++=a0=ifalt(*rl) /* templ=abs(Il) */ 

rle=Ql 

aO=-*rl 


*r2=a0=ifalt(*rl) 

/* temp2=abs(Ql) */ 


call _div (rl4) 
nop 

/* temp3=Q/I */ 


int24 temp2,templ,temp3 

.align 4 


call _atan (rl4) 
nop 

int24 tempi 
int24 temp3, theta 

/* find principal angle */ 

.align 4 

/* find arg(Z) where Z=Il+jQl */ 
rle=Il 



aO=*rl 
r2e=Q 1 

/* set flag if 1 1 is negative */ 


al=*r2 

r3e=theta 

r4e=pi 

/* set flag if Q1 is positive */ 


if(alt) pcgoto lefthalf 
nop 

/* 11 was negative! */ 


if(age) pcgoto quaddone 
nop 

/* Q1 positive - quadrant I */ 

quad4: 

*r3=aO=-*r3 
goto quaddone 
nop 

/* negate theta */ 

lefthalf: 

if(age) pcgoto quad2 
nop 

/* Q1 positive - quadrant II */ 

quad3: 

*r3=aO=*r3-*r4 
goto quaddone 
nop 

/* arg = theta - 180*/ 

quad2: 

*r3=a0=*r4-*r3 

/* arg = 1 80 - theta */ 


nop 

nop 

quaddone: rle = VpCoeff 

r2e = Vp 

*t2 = aO = *r3 * *rl /* Vp = VpCoeff * theta */ 

/****************************+********+*++i,++i l i t:ltt i ll # ilfltt *++++ it .**j. *******/ 

/* Variance = VarCoeff * lnlRl/R2l */ 

call _div (rl4) 
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nop 

int24 Rl,R2,templ 

.align 4 

call _loge (rl4) 
nop 

int24 temp3 
int24 templ,temp2 

.align 4 

rle = VarCoeff 

r2e = Var 

*r2 = aO = aO * *rl 

/*******************************************************************/ 
/* SNR = [IRll A (4/3)] / [RO * IR2l A (l/3) - IRll A (4/3)] */ 

/*******************************************************************y 

call _xtoy (rl4) /* R2 A (l/3) */ 

nop 

int24 temp3 

int24 R2,OneThird, tempi 

.align 4 

r8e = RO 

r9e = temp 1 /* *r9 = R2 A ( 1/3) */ 

rlOe = temp2 /* *rlO = Rl A (4/3) */ 

*r9 = aO = *r8 * aO /* *r9 = RO * R2 A (l/3) */ 

call _xtoy (rl4) 

nop 

int24 temp3 

int24 Rl,FourThirds,temp2 

.align 4 

*r9 = aO = *r9 - aO /* RO * R2 A ( 1/3) - Rl A (4/3) */ 

call _div (rl4) 

nop 

int24 temp2, tempi, SNR 

.align 4 

goto Out 
nop 

/* 

* Dump results to the output FIFO... 

*/ 

rsect ".Rl" 

Out: rle = DESTA 

r2e = QDATA 

r5 = 0x44 /* send results to board 4, zone 4 */ 

*rl = r5 

r5e = Vp /* mean doppler velocity */ 
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r7e = pack 

/* pack Vp and variance */ 

*r7++ = aO = int(*r5++) 

/* pack Vp */ 

*r7— = aO = int(*r5++) 

/* with variance...*/ 

r8e = Dwell_cnt 


r6 = *r8 

/* Get current dwell # */ 

nop 


aO = (*r2=*r7) + aO 

/* send packed result on Bus A */ 

r6 = r6 & 0x3f 

/* truncate to 6 bits */ 

r6 = r6»> 1 

/* shift right 1 bit */ 

r6 = r6»>l 

/* shift right 2 bits */ 

r6 = r6»> 1 

/* shift right 3 bits */ 

r6 = r6»> 1 

/* shift right 4 bits */ 

r6 = r6»> 1 

/* shift right 5 bits */ 

r6 = r6»>l 

/* shift right 6 bits */ 

r6 = r6»> 1 

/* shift right 7 bits */ 

r8e = Range_cnt 


00 

V-. 

* 

II 

/* get previous range # */ 

nop 


rl = rl + 1 

/* update it */ 

*r8 = rl 

/* and save it */ 

rl =rl &0x3ff 

/* truncate to 10 bits */ 

r6 = r6 1 rl 

/* and pack it into place */ 

*r7++ = aO = int(*r5) 

/* pack SNR */ 

nop 


*r7— = r6 

/* with marker */ 

nop 


nop 


nop 


aO = (*r2=*r7) + aO 

/* send SNR and marker on Bus A 

nop 


nop 


nop 


rle = RStart 


if(rl7— >=0) goto rl 

/* next range */ 

nop 


goto PassToken 

/* pass token to next process... */ 

nop 



/* 

* Pass the token A back to zone 0 
*/ 

PassToken: rle = DESTA 

r2e = QDATA 

r3 = Ox If /* pass token... */ 
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/* ...brd 1, zoned */ 


*rl = r3 
r3 = Ox 1 1 
*t2 = r3 

goto Dwell /* and prepare for next dwell */ 

nop 

.rsect ",R2" 

/ ***** * ***** * * * ****** * * * * % * ** * ********** * * * *************** 

/* Routine source file: atan.asm */ 

/* DSP32/DSP32C Application Software Library. Version 2.0 */ 

/ ********************** * * * * * * ************** * * * ** * * * * * ****** j 
/* 

_atan(A,B,C) 
float *A,*B,*C; 

scratch register short r4,r5 (for DSP32 only); 
scratch pointer register short rl,r2,r3,rl4; 
scratch register float a0,al,a2,a3; 

*/ 

#include <dspregs.h> /* Translates DSP32C keywords to DSP32 * 

* keywords. See Section 2.3.1 */ 

_atan: rle=*rl4++ /* pointer to local variable */ 

r2e=*rl4++ /* pointer to x */ 

r3e=_atanB 

*rl++=a3=*r2 + *r3 /* x+1 */ 

#if DSP32C 

a2 = seed(a3) 

*rl=a0=*r2 - *r3++ /* x-1 */ 

r2e = *rl4 /* pointer to out */ 

a0=*r3++ -a3 * a2 /* a0=1.4074-x*y */ 
al=*r3++ * a2 /*al=.81*x*/ 

a2=*r3++ -a3 * a2 /* a2=2.27424-x*y */ 

#e Ise 

*rl— =a0=*r2 - *r3++ /* x-1 */ 

r2e = *rl4 /* pointer to out */ 

nop 

r4=*rl++ /* load x+l=y into cau */ 

r5=*rl++ 

r4=r4 A 0xffff /* invert all bits except sign bit */ 

r5=r5 A Ox7fff 

*r2++=r4 /* 1st estimate of 1/y, 

or x, written to memory */ 

*r2-=r5 

a0=*r3++ -a3 * *r2 /* a0=1.4074-x*y */ 
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al=*r3++ * *r2 
a2=*r3++ -a3 * *r2 
#endif 

aO=*r3++ +aO * aO 

al=aO*al 

a2=a2*al 

nop 

a0=a2+a0*al 

nop 

nop 

al=*r3 + a3*a0 
aO = aO * *rl 
nop 
nop 

*r2 = a2 = aO-al*aO 

nop 

nop 

aO=a2*a2 

rle=_atanC 

nop 

a0=a0**r2 


/* al=.81*x */ 
f* a2=2.27424-x*y */ - 

/* a0= 1.71 742-2.8 148*x*y+x A 2*y A 2 */ 
/* al=1.14*x-.81*x A 2*y */ 

/* a2=1.842*x-.81*x A 2*y */ 


/* a0=3 . 8x-5 .4 1 x A 2 *y+3.42x A 3 * ... 
y A 2-.81x A 4*y A 3 */ 


/* al=-l + y* x(new) */ 
/* scale by "(x-1)" */ 


/* third iteration, ... */ 

/* a2=(x-l)/(x+l) */ 

/* y**2 */ 

/* y**3 */ 


al = *rl++ + aO * *rl++ /* c7 + c9*(y A 2) */ 

a2 = *rl++ + aO * *rl++ /* c3 + c5*(y A 2) */ 

aO = aO * *r2 /* y A 4 */ 

al = al * aO /* c7*(y A 3) + c9*(y A 5) */ 

a2 = *rl++ + a2 * aO /* pi/4 + c3*(y A 3) + c5*(y A 5) */ 


a2 = a2 + *rl++ * *r2 /* pi/4+cl*y+ ... 

c3*(y A 3)+c5*(y A 5) */ 

*r2 = aO = a2 + al * aO /* pi/4+cl*y+c3*(y A 3)+c5*(y A 5) 

... +c7*(y A 7)+c9*(y A 9) */ 

nop 

goto rl4+4 
nop 


_atanB: float 1 .0 

float 1.4074074347, 0.81, 2.27424702, -.263374728, -1.0 
_atanC: float .020835 1 , -.085 1 33 /* c9, c7 */ 

float. 180141, -.3302995 /* c5, c3 */ 

float .785398163, .999866 /* pi/4, cl */ 

/* */ 
/* _div(A,B,C) */ 
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/* */ 

/* float *A,*B,*C; - */ 

/* scratch register short r5,r6 (DSP32 only); */ 

/* scratch pointer register short rl,r2,r3,r4,rl4; */ 

/* scratch register float a0,al,a2,a3; */ 

/* */ 


/*******************************************************************/ 

_div: r4e = *rl4++ /* numerator pointer */ 

rle = *rl4++ /* denominator pointer */ 

r2e = *rl4 /* point to inverse location */ 

a3=*rl /* load number to be inverted into a3 */ 

#if DSP32C 

al = seed(*rl) 
r3e=_divA 
nop 

aO = *r3++ + a3 * al /* a0=- 1.601 4-x*y */ 
al = al * *r3++ /* al=.8237*x */ 

#else 

r3e=_divA 

r6=*rl++ /* load number into cau */ 

r5=*rl— 

r6=r6 A 0xffff /* invert all bits except sign bit */ 

r5=r5 A Ox7fff 

*r2++=r6 /* 1st estimate written to memory */ 

*r2— =r5 

aO = *r3++ + a3 * *r2 /* a0=- 1 ,6014-x*y */ 

al = *r2 * *r3++ /* al=.8237*x */ 

#endif 

nop 

aO = *r3++ + aO * aO /* a0=3.4166-3.2023*x* 

y+x A 2*y A 2 */ 

a2 = al * a3 /* a2=.8237*x*y */ 

nop 

aO = aO * al /* a0=2.814-2.638x A 2* 

y+.8237x A 3*y A 2 */ 

al = *r3 + aO * a2 /* al=1.814*x*y-2.638*x A 2* 

y A 2+.8237*x A 3*y A 3 */ 
nop 

aO = aO * *r4 /* scale by numerator */ 

nop 

nop 

*r2 = aO = aO - al * aO /* last iteration 

written to memory */ 
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nop 

goto rl4 + 4 
nop 


_divA: float -1.60138661 , 0.82371354 , 0.85221935 , -1.0 
.align 4 

/**********************************************************/ 
/* Routine source file: loge.asm */ 

/* DSP32/DSP32C Application Software Library. Version 2.0 */ 

/**********************************************************y 


/* 

_loge(A,B,C) 
float *A, *B, *C; 
scratch register short r5,r6; 
scratch pointer register short rl,r2,r3,r4,rl4; 
scratch register float a0,al,a2,a3; 

*/ 

#include <dspregs.h> /* Translates DSP32C keywords to DSP32 * 

* keywords. See Section 2.3.1 */ 


loge: r2e=*rl4++ 
rle=*rl4++ 
r4e=_logeC+36 
*r2=a0=*rl 
nop 

r3e=r2 + 4 


/* pointer to a copy of x */ 

/* pointer to operand, x */ 

/* pointer to c8 */ 

/* load x into copy */ 

/* pointer to scratch pad for exponent */ 


/* extract mantissa, M, and divide by 2 */ 
r5=-l + 128 /* -1 plus bias */ 

*r2=r51 /* replace exponent of operand with "-1" */ 

a3=*r2 /* move M/2 into a3 */ 


/* store exponent, E, as a floating point number */ 


r61=*rl 
rle=*rl4 
r6=r6+ 1-128 
*r3=r6 

*r3=aO=float(*r3) 


/* get exponent (byte) */ 

/* pointer to output */ 

/* remove bias and add 1 */ 

/* store as an integer, i.e., (E+l) */ 
/* convert to float and store */ 


/* compute log(M/2) in base 2 */ 
aO = a3 * a3 /* (M/2) A 2 */ 

al = *r4— + a3 * *r4— /*c7 + c8*M/2V 

a2 = *r4— + a3 * *r4- /* c4 + c5*M/2 */ 

aO = aO * a3 /* (M/2) A 3 */ 
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al = *r4— + al * a3 /* c6+c7*(M/2)+c8*((M/2) A 2) */ 

a3 = *r4— + a3 * *r4- /* cl + c2*M/2 */ 

a2 = *r4— + a2 * *r2 /* c3+(M/2)*(c4+c5*(M/2)) */ 

al = al * aO /* c6*((M/2) A 3)+:..+c8*((M/2) A 5) */ 

a3 = *r4-- + a3 * *r2 /* cO + c 1 *(M/2) + c2*((M/2) A 2) */ 

aO = *r3 * *r4 /* c9 * (E+l) */ 

a2 = a3 + a2 * aO /* cO+c 1 *(M/2)+...+c5*((M/2) A 5) */ 

al = a2 + al * aO /* cO+c 1 *(M/2)+...+c8*((M/2) A 8) */ 

/* add c9*(E +1) and store */ 

*rl = aO = al + aO 
nop 

goto rl4+4 
nop 

JogeC: float 0.693 1 47 1 806 /* c9 */ 
float -3.067466148 /* cO */ 

float 51.49518505 /* c3 */ 

float 11.30516183 /* cl */ 

float -27.7466647 /* c2 */ 

float -33.20167437 /* c6 */ 

float -66.69583126 /* c4 */ 

float 58.53503341 /* c5 */ 

float 10.98927014 /* c7 ♦/ 

float -1.613006222 /* c8 */ 

.align 4 


/* */ 

/* _sqr(A,B,C) */ 

/* */ 

/* float *A,*B,*C; */ 

/* scratch register short r4,r5; */ 

/* scratch pointer register short rl,r2,r3,r 14; */ 

/* scratch register float aO, a l,a2; */ 

/* */ 

/************* ************ ***********************************^*^^>^^^1 

#include <dspregs.h> /* Translates DSP32C keywords to DSP32 * 

* keywords. See Section 2.3.1 */ 


r2e = *rl4++ /* address of local variables */ 

rle = _sqrC 

*r2++ = aO = *rl++ /* stores exponent adjustment */ 

*r2— = aO = *rl /* stores exponent adjustment */ 


119 



rle = *rl4++ 
r3e = _sqrB 
r41 = *rl 
nop 

r5 = -r4 
*r2 = r51 
a2 = *rl * *r2 
r4 = r4 »1 
if (cc) pcgoto _sqrA 
al = a2 * a2 
r2e = r2 + 4 

_sqrA: aO = a2 * *r3++ 
aO = aO + *r3++ 
aO = aO + al * *r3++ 
al = a2 * *r3++ 
al = al + *r3++ 
aO = al + al * aO 
a2 = a2 * *r3++ 
r4 = r4 + 64 
aO = a2 * aO 
al = aO * aO 
*r2 = r41 
aO = aO * *r2 
al = *r3 - al * a2 
rle = *rl4 
nop 

*rl = aO = aO * al 
nop 

goto r 1 4+4 
nop 

_sqrB: float -0.38289166, 1.18787772,0.049693281,-1.92155474 
float 2.0667909, 0.5, 1.5000001 

_sqrC: float 1.0, 1.4142136 
.align 4 

/**********************************************************/ 
/* Routine source file: xtoy.asm */ 

/* DSP32/DSP32C Application Software Library. Version 2.0 */ 

/******************************************* ***************^ 

/* 

To use this routine with DSP32C, bit 4 of DAUC must be set to 0. 
_xtoy(A,B,C,D) 


/* address of input value */ 

/* load exponent into cau */ 

/* 2's complement inverts exponent */ 

/* reset exponent of operand to zero */ 

/* skip sqrt(2) adjustment */ 

/* x A 2 */ 

/* set pointer for sqrt(2) adjustment */ 

/* end of polynomial approximation */ 
/* start of iteration approximation */ 

/* final exponent = (opexp)/2 + 64 */ 

/* address of output location */ 

/* stores result */ 
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float *A, *B, *C, *D; 
scratch register short r5, r6; 
scratch pointer register short rl,r2,r3„r4,r7,rl4; 
scratch register float aO,al ,a2,a3; 

*/ 


#include <dspregs.h> 

/* Translates DSP32C keywords to DSP32 
* keywords. See Section 2.3.1 */ 

_xtoy: 

/* computes log(x) [base 2] */ 

r2e=*rl4++ 

/* pointer to a copy of x */ 

rle=*rl4++ 

/* pointer to operand, x */ 

r4e=_xtoyC+32 

/* pointer to c8 */ 

*r2=a0=*rl 

/* load x into copy */ 

r7e=*rl4++ 

/* pointer to y */ 

r3e=r2 + 4 

/* pointer to scratch pad for exponent */ 

/* extract mantissa, M, and divide by 2 */ 

r5=- 1 + 128 

/* -1 plus bias */ 

*r2=r51 

/* replace exponent of operand with "-1" */ 

a3=*r2 

/* move M/2 into a3 */ 

/* store exponent, E, as a 

floating point number */ 

r61=*rl 

nop 

/* get exponent (byte) */ 

r6=r6+l-128 

/* remove bias and add l */ 

*r3=r6 

/* store as an integer */ 

*r3=a0=float(*r3) 

/* convert to float and store */ 

/* compute log(M/2) in base 2 */ 

aO = a3 * a3 

/* (M/2) A 2 */ 

al = *r4— + a3 * *r4— 

/* c7 + c8*M/2 */ 

a2 = *r4— + a3 * *r4~ 

/* c4 + c5*M/2 */ 

aO = aO * a3 

/* (M/2) A 3 */ 

al = *r4— + al * a3 

/* c6 + c7*(M/2) + c8*((M/2) A 2) */ 

a3 = *r4 — H a3 * *r4— 

/* cl + c2*M/2 */ 

a2 = *r4— + a2 * *r2 

/* c3 + (M/2)*(c4 + c5*(M/2)) */ 

al = al * aO 

/* c6*((M/2) A 3)+c7*((M/2) A 4)+ 
c8*((M/2) A 5) */ 

a3 = *r4 + a3 * *r2 

/* cO + cl*(M/2) + c2*((M/2) A 2) */ 

a2 = a3 + a2 *a0 

/* cO + cl*(M/2) + ... + 
c5*((M/2) A 5) */ 

al = a2 + al * aO 

/* cO + cl*(M/2) + ... + 
c8*((M/2) A 8) */ 

aO = al + *r3 

I* add (E+l) ==> aO=log(x) [base 2] *! 


/* computes 2 A [y*log(x)} [base 2] */ 
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r2e = *rl4++ 

/* 

r3e = _xtoyD 

/* 

aO = aO * *r7 

/* 


r2 points to output location */ 
r3 points to 0.5 and then ci’s */ 
y*log(x) [base 2] */ 


a3 = a0 - *r3++ /* x - 0.5 */ 

*r2 = a3 = int(a3) /* stores i, the "integer" part */ 


a3 = float(a3) 
a3 = -a3 + aO 


/* converts i to float format */ 

/* computes f, the "fractional part */ 


nop 

r61 = *r2 


/* stores i in r61 */ 


/* computes 2 A f below */ 
aO = a3 * a3 

al = *r3++ + a3 * *r3++ 
a2 = *r3++ + a3 * *r3++ 
aO = aO * aO 
a3 = *r3 + a3 * *r3++ 
al = a3 + al * aO 
*r2 = a2 = al + a2 * aO 

r6 = r6+ 128 

nop 

nop 

*r2 = r61 
aO = *r2 
return (r 14) 
nop 


/* f A 2 */ 

/* c2+c3(f) */ 

/* c4+c5(f) */ 

/* f A 4 */ 

/* cO+cl(f) */ 

/* cO+c 1 (f)+c2(f A 2)+c3(f A 3) */ 
/* cO+cl(f)+...+c5(f A 5) */ 

/* adjusts exponent for bias */ 
/* waits for DAU to */ 

/* write the mantissa */ 

/* merge in exponent */ 


xtoyC: float -4.4254182 

/* 

cO 

*/ 

float 74.29184809343 

/* 

c3 

*/ 

float 16.3099009156 

/* 

cl 

*/ 

float -40.02997556826 

/* 

c2 

*/ 

float -47.89989096077 

/* 

c6 

*/ 

float -96.221745 

/* 

c4 

*/ 

float 84.44820241827 

/* 

c5 

*/ 

float 15.8541655496 

/* 

c7 

*/ 

float -2.32707607725 

/* 

c8 

*/ 


_xtoyD: float 0.5 

float 0.0558263 1 80623292 /* c3 */ 

float 0.240153617040129 /* c2 */ 

float 0.00 1 8775766770276 /* c5 */ 

float 0.00898934008333 12 /* c4 */ 
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float 0.6931530732007278 /* cl */ 

float 0.9999999702 /* c 0 */ 


.rsect ",R0" 
.align 4 

IpPrev: 

2*float 42.0 

/* Contains Ip from previous block */ 

Ip: 

200*float 

/* Contains Ip for current block */ 


float 0.0, 0.0 

/* #pulses each range */ 

/* Overshoot zone. DON'T REMOVE! */ 

.rsect M .R1" 
.align 4 

QpPrev: 

2*float 69.0 

/* Contains Qp from previous block */ 

Qp: 

200*float 

/* Contains Qp for current block */ 


float 0. 0,0.0 

/* #pulses each range */ 

/* Overshoot zone. DON'T REMOVE! */ 

VpCoeff: 

float 5.2521 131 

/* Coefficient needed to get mean Doppler */ 

VarCoeff: 

float 18.389794 

/* VpCoeff = (lambda/2)*[PRF/(2*pi)] */ 
/* Ditto for variance of the time series */ 

OneThird: 

float 1. 0/3.0 

/* VarCoeff = lambda A 2*PRF A 2/[24*pi A 2] */ 

FourThirds: 

float 4.0/3.0 


pi: 

float 3.1415926535898 

Zero: 

float 0.0 


temp 1 : 

float 0.0 

/* scratch registers... */ 

temp2: 

float 0.0 

temp3: 

float 0.0 


temp4: 

float 0.0 


Routaddr: 

int24 0 

/* scratch address register */ 

.align 4 

Count 1: 

int 598 

/* #range gates * #pulses each range -2*1 

Count2: 

int 199 

/* Note: #pulses means ONE BLOCK ONLY! ! 
/* #pulses each range - 1 */ 

Count3: 

int 10 

/* again, #pulses ONE BLOCK ONLY! */ 
/* Distance between consequtive pulses in 

Count4: 

int 26 

* packed I & Q array for a given range 

* distance = #range gates * 4 - 2 */ 

/* (#range gates * 9) - 1 */ 

NumBlock: 

int 8 

/* # blocks of data - 2 that comprise a Dwell */ 

RangeGates: 

int 1 

/* Number of range gates -2*1 

.align 4 

NumPulse: 

float 1998.0 

/* Number of Pulses each range - 2 


* INCLUDES ALL BLOCKS (IE. TOTAL 

* # PULSES - 2 FOR ALL BLOCKS !!!)*/ 
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pack: 

float 0.0 

/* 

holdpack: 

int 0 

/* 

PackFlag: 

int 0 

/* 

* 

Fubar: 

int Oxaaaa 

* 

/* 

stopflag: 

int 1 

/* 

addrflag: 

float 0.0 

/* 

.rsect ".R2" 



.align 4 



10: 

float 0.0 

/* 

11: 

float 0.0 

/* 

12: 

float 0.0 

/* 

Ql: 

float 0.0 

/* 

Q2: 

float 0.0 


theta: 

float 0.0 

/* 

R0: 

float 0.0 

/* 

Rl: 

float 0.0 

/* 

R2: 

float 0.0 

/* 

Vp: 

float 0.0 

/* 

Var: 

float 0.0 

/* 

SNR: 

float 0.0 

/* 


.rsect ".text" 
.align 4 
IQPack: 

1200*int Oxffff 

/* 



/* 

.align 4 
PartSum: 

27*float 23.0 

/* 



/* 

.align 4 
Range_cnt: 

int 0 

/* 

Dwell cnt: 

int 0 

/* 


holding spot for result packing */ 
holds extra result between range packing */ 
>0 => SNR from last range waiting to be 
packed. It is sitting in holdpack. 

=0 => previous range completely packed */ 
For test purposes only */ 

For test purposes only */ 
for test purposes only */ 


IMPORTANT!! The variables I0-Q2 */ 
*MUST* remain in this order - the */ 
program assumes they are blocked as */ 
shown when loading them... */ 

argument of lag 1 autocorrelation */ 
magnitude of lag 0 autocorrelation */ 
magnitude of lag 1 autocorrelation */ 
magnitude of lag 2 autocorrelation */ 
Doppler velocity */ 

Variance */ 

Signal to Noise Ratio */ 


buffer holding packed I & Q */ 
2*#range gates*#pulses each range */ 

Partial I & Q sums */ 
should be 9*#range gates */ 

for test purposes */ 
for test purposes */ 
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